Add the TRAP-chip simulation by Clemens
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 4 Apr 2008 13:00:39 +0000 (13:00 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 4 Apr 2008 13:00:39 +0000 (13:00 +0000)
TRD/AliTRDfeeParam.h
TRD/AliTRDmcmSim.cxx
TRD/AliTRDmcmSim.h
TRD/AliTRDrawData.cxx
TRD/AliTRDrawData.h
TRD/AliTRDtrapAlu.cxx [new file with mode: 0644]
TRD/AliTRDtrapAlu.h [new file with mode: 0644]

index 50cd1a5..3fe55c6 100644 (file)
@@ -71,6 +71,17 @@ class AliTRDfeeParam : public TObject
   static  Int_t    GetPFtimeConstant()    { return fgkPFtimeConstant;   }
   static  Int_t    GetPFeffectPedestal()  { return fgkPFeffectPedestal; }
 
+  //new
+  static  Int_t    GetQacc0Start()        {  return fgkPREPqAcc0Start; }
+  static  Int_t    GetQacc0End()          {  return fgkPREPqAcc0End; }
+  static  Int_t    GetQacc1Start()        {  return fgkPREPqAcc1Start; }
+  static  Int_t    GetQacc1End()          {  return fgkPREPqAcc1End; }
+          Float_t  GetMinClusterCharge() const {  return fgkMinClusterCharge; }
+  static  Int_t    GetLinearFitStart()    {  return fgkPREPLinearFitStart; }
+  static  Int_t    GetLinearFitEnd()      {  return fgkPREPLinearFitEnd;  }
+
+
+
   //        Float_t  GetClusThr()           { return fClusThr; };
   //        Float_t  GetPadThr() const { return fPadThr; };
   //        Int_t    GetTailCancelation() const { return fTCOn; };
@@ -85,6 +96,10 @@ class AliTRDfeeParam : public TObject
           Float_t  GetTFc1()        const { return fTFc1;           }
           Float_t  GetTFc2()        const { return fTFc2;           }
 
+ // for tracklets
+         Bool_t   GetTracklet()         const { return fgkTracklet; } 
+          Int_t    GetMaxNrOfTracklets() const { return fgkMaxNrOfTracklets; } 
+
   static  Float_t  GetTFattPar()          { return ((Float_t) fgkTFattPar1) / ((Float_t) fgkTFattPar2); }
           Float_t  GetTFf0()        const { return 1.0 + fgkTFon*(-1.0+GetTFattPar()); }   // 1 if TC off
 
@@ -145,11 +160,17 @@ class AliTRDfeeParam : public TObject
   static const Int_t    fgkGFnoise           =  0;        // Noise level increased by gain filter x 100 [in ADC] (to be measured)
 
   // TF setup
-  static const Int_t    fgkTFtype            = 1;         // TC type (0=analog, 1=digital, 2=MI) (old name fFilterType)
+  static const Int_t    fgkTFtype            = 1;         // TC type (0=analog, 1=digital, 2=MI, 3=close to electronics) (old name fFilterType)
 
   // OLD TF setup (calculated from above)  (valid only for fgkTFsimType = 0 or 1)
   static const Int_t    fgkTFnExp          = 1;           // Number of exponential for simType 0 and 1
 
+ // Tracklet  processing on/off 
+  static const Bool_t   fgkTracklet         = kFALSE; // default should be kFALSE
+
+  // max. nr of tracklet words for one mcm
+  static const Int_t    fgkMaxNrOfTracklets = 4; 
+
   // following need Instance because initialized in constructor
                Float_t  fTFr1;                            // Time constant [us] long (old name fR1)
                Float_t  fTFr2;                            // Time constant [us] short(old name fR2)
@@ -170,12 +191,19 @@ class AliTRDfeeParam : public TObject
                Int_t    fEBignoreNeighbour;               // EBIN 0:include neighbor
 
   // Charge accumulators
-  static const Int_t    fgkPREPqAcc0Start     =  0;       // Preprocessor Charge Accumulator 0 Start
+  static const Int_t    fgkPREPqAcc0Start     =  5;       // Preprocessor Charge Accumulator 0 Start
   static const Int_t    fgkPREPqAcc0End       = 10;       // Preprocessor Charge Accumulator 0 End
   static const Int_t    fgkPREPqAcc1Start     = 11;       // Preprocessor Charge Accumulator 1 Start
   static const Int_t    fgkPREPqAcc1End       = 20;       // Preprocessor Charge Accumulator 1 End
   static const Int_t    fgkMinClusterCharge   = 20;       // Hit detection [in ADC units]
 
+  //new
+  //time constants for linear fit
+  static const Int_t    fgkPREPLinearFitStart = 5;
+  static const Int_t    fgkPREPLinearFitEnd   = 20;
+
+
+
   // OLD TRAP processing parameters calculated from above
   //static const Float_t  fClusThr;                       // Cluster threshold
   //static const Float_t  fPadThr;                        // Pad threshold
index b044b66..292235c 100644 (file)
@@ -80,6 +80,16 @@ The default raw version is 2.
                                                                  Ken Oyama
 */
 
+// if no histo is drawn, these are obsolete
+#include <TH1.h>
+#include <TCanvas.h>
+
+// only needed if I/O of tracklets is activated
+#include <TObject.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TSystem.h>
+
 #include <fstream>
 
 #include <TMath.h>
@@ -92,11 +102,18 @@ The default raw version is 2.
 #include "AliTRDgeometry.h"
 #include "AliTRDcalibDB.h"
 
+// additional for new tail filter and/or tracklet
+#include "AliTRDtrapAlu.h"
+#include "AliTRDpadPlane.h"
+
+
 ClassImp(AliTRDmcmSim)
 
 //_____________________________________________________________________________
 AliTRDmcmSim::AliTRDmcmSim() :TObject()
   ,fInitialized(kFALSE)
+  ,fNextEvent(-1)    //new
+  ,fMaxTracklets(-1) //new
   ,fChaId(-1)
   ,fSector(-1)
   ,fStack(-1)
@@ -108,6 +125,9 @@ AliTRDmcmSim::AliTRDmcmSim() :TObject()
   ,fRow(-1)
   ,fADCR(NULL)
   ,fADCF(NULL)
+  ,fADCT(NULL)      //new
+  ,fPosLUT(NULL)    //new
+  ,fMCMT(NULL)      //new
   ,fZSM(NULL)
   ,fZSM1Dim(NULL)
   ,fFeeParam(NULL)
@@ -126,7 +146,9 @@ AliTRDmcmSim::AliTRDmcmSim() :TObject()
 //_____________________________________________________________________________
 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m) 
   :TObject(m)
-  ,fInitialized(kFALSE)
+  ,fInitialized(kFALSE) 
+  ,fNextEvent(-1)    //new
+  ,fMaxTracklets(-1) //new
   ,fChaId(-1)
   ,fSector(-1)
   ,fStack(-1)
@@ -138,12 +160,16 @@ AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
   ,fRow(-1)
   ,fADCR(NULL)
   ,fADCF(NULL)
+  ,fADCT(NULL)      //new
+  ,fPosLUT(NULL)    //new
+  ,fMCMT(NULL)      //new
   ,fZSM(NULL)
   ,fZSM1Dim(NULL)
   ,fFeeParam(NULL)
   ,fSimParam(NULL)
   ,fCal(NULL)
   ,fGeo(NULL)
+  
 {
   //
   // AliTRDmcmSim copy constructor
@@ -164,13 +190,21 @@ AliTRDmcmSim::~AliTRDmcmSim()
     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
       delete [] fADCR[iadc];
       delete [] fADCF[iadc];
+      delete [] fADCT[iadc];
       delete [] fZSM [iadc];
     }
     delete [] fADCR;
     delete [] fADCF;
+    delete [] fADCT;
     delete [] fZSM;
     delete [] fZSM1Dim;
   }
+  if(fInitialized){
+    delete [] fPosLUT;
+    delete [] fMCMT;
+  }
+
   delete fGeo;
 
 }
@@ -195,7 +229,8 @@ void AliTRDmcmSim::Copy(TObject &m) const
   //
   // Copy function
   //
-
+  ((AliTRDmcmSim &) m).fNextEvent     = 0; //new
+  ((AliTRDmcmSim &) m).fMaxTracklets  = 0; //new
   ((AliTRDmcmSim &) m).fInitialized   = 0;
   ((AliTRDmcmSim &) m).fChaId         = 0;
   ((AliTRDmcmSim &) m).fSector        = 0;
@@ -208,6 +243,9 @@ void AliTRDmcmSim::Copy(TObject &m) const
   ((AliTRDmcmSim &) m).fRow           = 0;
   ((AliTRDmcmSim &) m).fADCR          = 0;
   ((AliTRDmcmSim &) m).fADCF          = 0;
+  ((AliTRDmcmSim &) m).fADCT          = 0; //new
+  ((AliTRDmcmSim &) m).fPosLUT        = 0; //new
+  ((AliTRDmcmSim &) m).fMCMT          = 0; //new
   ((AliTRDmcmSim &) m).fZSM           = 0;
   ((AliTRDmcmSim &) m).fZSM1Dim       = 0;
   ((AliTRDmcmSim &) m).fFeeParam      = 0;
@@ -218,13 +256,15 @@ void AliTRDmcmSim::Copy(TObject &m) const
 }
 
 //_____________________________________________________________________________
-void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
+//void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent ) // only for readout tree (new event)
+void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos ) // only for readout tree (new event)
 {
   //
   // Initialize the class with new geometry information
   // fADC array will be reused with filled by zero
   //
 
+  fNextEvent     = 0; //**!!new!!**
   fFeeParam      = AliTRDfeeParam::Instance();
   fSimParam      = AliTRDSimParam::Instance();
   fCal           = AliTRDcalibDB::Instance();
@@ -239,15 +279,28 @@ void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
   fNTimeBin      = fCal->GetNumberOfTimeBins();
   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
 
+  fMaxTracklets  = fFeeParam->GetMaxNrOfTracklets();
+
+
+  //if ((fChaId == 0 && fRobPos == 0 && fMcmPos == 0)) {
+  /* if (newEvent == kTRUE) {
+    fNextEvent = 1;
+    }*/
+   fNextEvent = 0;
+
+
   // Allocate ADC data memory if not yet done
   if( fADCR == NULL ) {
     fADCR    = new Int_t *[fNADC];
     fADCF    = new Int_t *[fNADC];
+    fADCT    = new Int_t *[fNADC]; //new
     fZSM     = new Int_t *[fNADC];
     fZSM1Dim = new Int_t  [fNADC];
     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
       fADCR[iadc] = new Int_t[fNTimeBin];
       fADCF[iadc] = new Int_t[fNTimeBin];
+      fADCT[iadc] = new Int_t[fNTimeBin]; //new
       fZSM [iadc] = new Int_t[fNTimeBin];
     }
   }
@@ -257,11 +310,24 @@ void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
       fADCR[iadc][it] = 0;
       fADCF[iadc][it] = 0;
+      fADCT[iadc][it] = -1;  //new
       fZSM [iadc][it] = 1;   // Default unread = 1
     }
     fZSM1Dim[iadc] = 1;      // Default unread = 1
   }
   
+  //new:
+  fPosLUT = new Int_t[128];
+  for(Int_t i = 0; i<128; i++){
+    fPosLUT[i] = 0;
+  }
+  
+  fMCMT = new UInt_t[fMaxTracklets];
+  for(Int_t i = 0; i < fMaxTracklets; i++) {
+    fMCMT[i] = 0;
+  }
+
+
   fInitialized = kTRUE;
 }
 
@@ -279,6 +345,125 @@ Bool_t AliTRDmcmSim::CheckInitialized()
 }
 
 //_____________________________________________________________________________
+
+
+void AliTRDmcmSim::SetPosLUT() {
+  Double_t iHi = (Double_t)fCal->GetPRFhi();
+  Double_t iLo = (Double_t)fCal->GetPRFlo();
+  Int_t   nBin = fCal->GetPRFbin();
+  Int_t   iOff = fLayer * nBin;
+  Int_t kNplan = fGeo->Nplan();
+
+  Float_t  *sPRFsmp   = new Float_t[nBin*kNplan];
+  Double_t *sPRFlayer = new Double_t[nBin];
+  
+  
+  for(Int_t i = 0; i<nBin*kNplan; i++){
+    
+    //printf("%f\n",fCal->GetSampledPRF()[i]);
+    sPRFsmp[i] = fCal->GetSampledPRF()[i]; 
+  
+  }
+
+  Double_t sWidth = (iHi-iLo)/((Double_t) nBin);
+  Int_t   sPad    = (Int_t) (1.0/sWidth);
+  
+  // get the PRF for actual layer (interpolated to ibin data-points; 61 measured)
+  for(Int_t iBin = 0; iBin < nBin; iBin++){
+    sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin];
+  }
+
+  Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5);                           // bin-nr. for pad-position 0
+  
+  Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5);          // bin-nr. for pad-position 0.5
+  bin1 = bin1 + 1;
+  bin0 = bin0 + 1;  //avoid negative values in aYest (start right of symmetry center)
+  while (bin0-sPad<0) {
+    bin0 = bin0 + 1;
+  }
+  while (bin1+sPad>=nBin) {
+    bin1 = bin1 - 1;
+  }
+  
+  Double_t* aYest = new Double_t[bin1-bin0+1];
+
+  /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5);
+  TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5);
+  TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5);
+  TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5);
+  TCanvas *c1 = new TCanvas("c1","c1",800,1000);
+  hist1->Draw();
+  TCanvas *c2 = new TCanvas("c2","c2",800,1000);
+  hist2->Draw();
+  TCanvas *c3 = new TCanvas("c3","c3",800,1000);
+  hist3->Draw();
+  TCanvas *c4 = new TCanvas("c4","c4",800,1000);
+  hist4->Draw();*/
+  
+  for(Int_t iBin = bin0; iBin <= bin1; iBin++){
+    aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1
+    //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo;
+    //  hist1->Fill(position,aYest[iBin-bin0]);
+  }
+  
+
+
+  Double_t aY[128]; // reversed function
+
+  AliTRDtrapAlu a;
+  a.Init(1,8,0,31);
+  
+  for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries; 
+    Double_t yest = ((Double_t)j)/256; 
+    
+    Int_t iBin = 0;
+    while (yest>aYest[iBin] && iBin<(bin1-bin0)) {
+      iBin = iBin+1;
+    }
+    if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) {
+      aY[j] = 0.5;                      // yest too big
+      //hist2->Fill(yest,aY[j]);
+      
+    }
+    else {
+      Int_t bin_d = iBin + bin0 - 1;
+      Int_t bin_u = iBin + bin0;
+      Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y
+      Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y
+      Double_t yest_d = aYest[iBin-1];                     // lower estimated y
+      Double_t yest_u = aYest[iBin];                       // upper estimated y
+      
+      aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d;
+      //hist2->Fill(yest,aY[j]);
+     
+    }
+    aY[j] = aY[j] - yest;
+    //hist3->Fill(yest,aY[j]);
+    // formatting
+    a.AssignDouble(aY[j]);
+    //a.WriteWord();
+    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
+    //hist4->Fill(yest,fPosLUT[j]);
+    
+  }
+   
+  
+  delete [] sPRFsmp;
+  delete [] sPRFlayer;
+  delete [] aYest;
+  
+}
+
+
+//_____________________________________________________________________________
+Int_t* AliTRDmcmSim::GetPosLUT(){
+  return fPosLUT;
+}
+
+
+
 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
 {
   //
@@ -420,6 +605,40 @@ Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
 }
 
 //_____________________________________________________________________________
+Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
+{
+  //
+  // Produce tracklet data stream from this MCM and put in buf
+  // Returns number of words filled, or negative value 
+  // with -1 * number of overflowed words
+  //
+
+  UInt_t  x;
+  Int_t   nw  = 0;  // Number of written words
+  Int_t   of  = 0;  // Number of overflowed words
+    
+  if( !CheckInitialized() ) return 0;
+
+  // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM 
+  // fMCMT is filled continuously until no more tracklet words available
+
+  Int_t wd = 0;
+  while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
+      x = fMCMT[wd];
+      if (nw < maxSize) {
+       buf[nw++] = x;
+      }
+      else {
+       of++;
+      }
+      wd++;
+  }
+  
+  if( of != 0 ) return -of; else return nw;
+}
+
+
+//_____________________________________________________________________________
 void AliTRDmcmSim::Filter()
 {
   //
@@ -510,6 +729,18 @@ void AliTRDmcmSim::FilterTail()
       }
     }
     break;
+
+    //new
+  case 3: // Exponential filter using AliTRDtrapAlu class
+    for (Int_t iCol = 0; iCol < fNADC; iCol++) {
+      FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
+      for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
+       fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
+       fADCT[iCol][iTime] = itarg[iTime];    // 12bits; to be used for tracklet; tracklet will have own container; 
+      }
+    }
+    break;
+
     
   default:
     AliError(Form("Invalid filter type %d ! \n", tftype ));
@@ -937,3 +1168,1570 @@ void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
   }
 }
 
+
+//_____________________________________________________________________________________
+//the following filter uses AliTRDtrapAlu-class
+
+void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
+  //static Int_t count = 0;
+  Double_t dt = 0.1;
+  Double_t r1 = (Double_t)fFeeParam->GetTFr1();
+  Double_t r2 = (Double_t)fFeeParam->GetTFr2();
+  Double_t c1 = (Double_t)fFeeParam->GetTFc1();
+  Double_t c2 = (Double_t)fFeeParam->GetTFc2();
+  
+  nexp = 1;
+
+  //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
+  //parameters need to be adjusted
+  AliTRDtrapAlu lambdaL;
+  AliTRDtrapAlu lambdaS;
+  AliTRDtrapAlu alphaL;
+  AliTRDtrapAlu alphaS;
+  
+  AliTRDtrapAlu correction;
+  AliTRDtrapAlu result;
+  AliTRDtrapAlu bufL;
+  AliTRDtrapAlu bufS;
+  AliTRDtrapAlu bSource;
+  
+  lambdaL.Init(1,11);
+  lambdaS.Init(1,11);
+  alphaL.Init(1,11);
+  alphaS.Init(1,11);
+  
+  //count=count+1;
+
+  lambdaL.AssignDouble(TMath::Exp(-dt/r1));
+  lambdaS.AssignDouble(TMath::Exp(-dt/r2));
+  alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
+  alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
+  
+  //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
+  correction.Init(10,2);
+  result.Init(10,2);
+  bufL.Init(10,2);
+  bufS.Init(10,2);
+  bSource.Init(10,2);
+  
+  for(Int_t i = 0; i < n; i++) {
+    bSource.AssignInt(source[i]);
+    result = bSource - correction; // subtraction can produce an underflow
+    if(result.GetSign() == kTRUE) {
+      result.AssignInt(0);
+    }
+    
+    //target[i] = result.GetValuePre();  // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
+    
+    target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2 
+
+    //printf("target-Wert zur Zeit %d : %d",i,target[i]);
+    //printf("\n");
+    
+    bufL  =  bufL + (result * alphaL);
+    bufL  =  bufL * lambdaL; 
+    
+    bufS  =  bufS + (result * alphaS);
+    bufS  =  bufS * lambdaS;  // eventually this should look like:
+    // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
+
+    correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
+  }
+  
+}
+
+
+
+
+
+
+
+
+//__________________________________________________________________________________
+
+
+// in order to use the Tracklets, please first 
+// -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
+// -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
+// currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
+
+// 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
+
+void AliTRDmcmSim::Tracklet(){
+    // tracklet calculation
+    // if you use this code outside a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
+
+  if(!CheckInitialized()){ return; }
+  
+  Bool_t filtered = kTRUE;
+  
+  
+  
+  AliTRDtrapAlu data;
+  data.Init(10,2);
+  if(fADCT[0][0]==-1){                      // check if filter was applied
+    filtered = kFALSE;
+    for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
+      for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
+       data.AssignInt(fADCR[iadc][iT]);
+       fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called this is done correctly
+      }
+    }
+   
+  }
+  
+  // the online ordering of mcm's is reverse to the TRAP-ordering(?)! reverse fADCT (to be consistent to TRAP), then do all calculations
+  // reverse fADCT:
+  Int_t** rev0 = new Int_t *[fNADC];
+  Int_t** rev1 = new Int_t *[fNADC];
+  
+  for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
+    rev0[iadc] = new Int_t[fNTimeBin];
+    rev1[iadc] = new Int_t[fNTimeBin];
+    for( Int_t iT = 0; iT < fNTimeBin; iT++) {
+      if( iadc <= fNADC-iadc-1 ) {
+       rev0[iadc][iT]  = fADCT[fNADC-iadc-1][iT];
+       rev1[iadc][iT]  = fADCT[iadc][iT];
+       fADCT[iadc][iT] = rev0[iadc][iT];
+      }
+      else {
+       rev0[iadc][iT]  = rev1[fNADC-iadc-1][iT];
+       fADCT[iadc][iT] = rev0[iadc][iT];
+      }
+    }
+  }
+  for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
+    delete[] rev0[iadc];
+    delete[] rev1[iadc];
+  }
+  
+  delete[] rev0;
+  delete[] rev1;
+  
+  rev0 = NULL;
+  rev1 = NULL;
+    
+  // get the filtered pedestal; supports only electronic tail-cancellation filter
+  AliTRDtrapAlu filPed;
+  Int_t ep = 0;
+  Int_t *ieffped = new Int_t[fNTimeBin];
+  for(Int_t iT = 0; iT < fNTimeBin; iT++){
+    ieffped[iT] = ep; 
+  }
+  if( filtered == kTRUE ) {
+    if( fFeeParam->IsPFon() ){
+      ep = fFeeParam->GetPFeffectPedestal();
+    }
+    Int_t      nexp  = fFeeParam->GetTFnExp();
+    Int_t  *isource  = new Int_t[fNTimeBin];
+    filPed.Init(10,2);
+    filPed.AssignInt(ep);           
+    Int_t epf = filPed.GetValue();  
+    for(Int_t iT = 0; iT < fNTimeBin; iT++){
+      isource[iT] = ep;                  
+      ieffped[iT] = epf;
+    }
+    if( fFeeParam->IsTFon() ) {
+      FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
+    }
+  
+    delete[] isource;
+  }
+  
+  //the following values should be in AliTRDfeeParam and have to be read in properly
+  //naming follows conventions in TRAP-manual
+  
+  
+  Bool_t bVBY = kTRUE;                         // cluster-verification bypass
+
+  Double_t cQTParam = 0;                      // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
+  AliTRDtrapAlu cQTAlu; 
+  cQTAlu.Init(1,10,0,63);
+  cQTAlu.AssignDouble(cQTParam);
+  Int_t cQT = cQTAlu.GetValue();
+
+  // linear fit 
+  Int_t tFS = fFeeParam->GetLinearFitStart();  // linear fit start
+  Int_t tFE = fFeeParam->GetLinearFitEnd();    // linear fit stop
+   
+  // charge accumulators
+  Int_t tQS0 = fFeeParam->GetQacc0Start();     // start-time for charge-accumulator 0
+  Int_t tQE0 = fFeeParam->GetQacc0End();       // stop-time for charge-accumulator 0
+  Int_t tQS1 = fFeeParam->GetQacc1Start();     // start-time for charge-accumulator 1 
+  Int_t tQE1 = fFeeParam->GetQacc1End();       // stop-time for charge-accumulator 1
+  // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
+  Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
+  AliTRDtrapAlu cTHAlu;  
+  cTHAlu.Init(12,2);
+  cTHAlu.AssignDouble(cTHParam);
+  Int_t cTH = cTHAlu.GetValue();                                 // cTH used for comparison
+
+  struct List_t {
+    List_t *next;
+    Int_t iadc;
+    Int_t value;
+  };
+  
+  List_t selection[7];            // list with 7 elements
+  List_t *list = NULL;
+  List_t *listLeft = NULL;
+    
+  Int_t* qsum = new Int_t[fNADC];
+   
+  // fit sums
+  AliTRDtrapAlu qsumAlu;
+  qsumAlu.Init(12,2);           // charge sum will be 12+2 bits
+  AliTRDtrapAlu dCOGAlu; 
+  dCOGAlu.Init(1,7,0,127);      // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
+  AliTRDtrapAlu yrawAlu;
+  yrawAlu.Init(1,8,-1,255);
+  AliTRDtrapAlu yAlu;
+  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 -
+  AliTRDtrapAlu xAlu;
+  xAlu.Init(5,8);               // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
+  AliTRDtrapAlu xxAlu;
+  xxAlu.Init(10,0);            
+  AliTRDtrapAlu yyAlu;
+  yyAlu.Init(1,16,0,0xFFFF);    // maximum is 2^16-1; 16Bit for past-commas
+  AliTRDtrapAlu xyAlu;
+  xyAlu.Init(6,8);
+  AliTRDtrapAlu XAlu;
+  XAlu.Init(9,0);
+  AliTRDtrapAlu XXAlu;
+  XXAlu.Init(14,0);
+  AliTRDtrapAlu YAlu;
+  YAlu.Init(5,8);               // 14 bit, 1 is sign-bit; therefore only 13 bit 
+  AliTRDtrapAlu YYAlu;
+  YYAlu.Init(5,16);
+  AliTRDtrapAlu XYAlu;
+  XYAlu.Init(8,8);              // 17 bit, 1 is sign-bit; therefore only 16 bit        
+  AliTRDtrapAlu qtruncAlu;
+  qtruncAlu.Init(12,0);
+  AliTRDtrapAlu QT0Alu;
+  QT0Alu.Init(15,0);
+  AliTRDtrapAlu QT1Alu;
+  QT1Alu.Init(16,0);
+  
+  AliTRDtrapAlu inverseNAlu;
+  inverseNAlu.Init(1,8);        // replaces the LUT for 1/N
+  AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
+  MeanChargeAlu.Init(8,0);
+  AliTRDtrapAlu TotalChargeAlu;
+  TotalChargeAlu.Init(17,8);
+  //nr of post comma bits should be the same for inverseN and TotalCharge
+  
+  
+  SetPosLUT();                    // initialize the position correction LUT for this MCM;
+
+
+  // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
+  Int_t *X   = new Int_t[fNADC-2];
+  Int_t *XX  = new Int_t[fNADC-2];
+  Int_t *Y   = new Int_t[fNADC-2];
+  Int_t *YY  = new Int_t[fNADC-2];
+  Int_t *XY  = new Int_t[fNADC-2];
+  Int_t *N   = new Int_t[fNADC-2];
+  Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
+  Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
+  
+  for (Int_t iCol = 0; iCol < fNADC-2; iCol++) { 
+      
+      // initialize fit-sums 
+      X[iCol]   = 0;
+      XX[iCol]  = 0;
+      Y[iCol]   = 0;
+      YY[iCol]  = 0;
+      XY[iCol]  = 0;
+      N[iCol]   = 0;
+      QT0[iCol] = 0;
+      QT1[iCol] = 0;
+  }
+  
+
+  filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
+  
+  for(Int_t iT = 0; iT < fNTimeBin; iT++){
+    
+    if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no? // !!**enable**!!
+
+    // reset
+    Int_t portChannel[4]   = {-1,-1,-1,-1};   
+    Int_t clusterCharge[4] = {0,0,0,0};
+    Int_t leftCharge[4]    = {0,0,0,0};
+    Int_t centerCharge[4]  = {0,0,0,0}; 
+    Int_t rightCharge[4]   = {0,0,0,0};
+    
+    Int_t mark = 0;
+    
+    filPed.AssignFormatted(ieffped[iT]);   // no size-checking with AssignFormatted; ieffped>=0
+    filPed = filPed;                       // this checks the size
+    
+    ieffped[iT] = filPed.GetValue();
+        
+    for(Int_t i = 0; i<7; i++){
+      selection[i].next       = NULL;
+      selection[i].iadc       =   -1;     // value of -1: invalid adc
+      selection[i].value      =    0;
+   
+    }
+    // selection[0] is starting list-element; just for pointing
+
+    // loop over inner adc's 
+    for (Int_t iCol = 1; iCol < fNADC-1; iCol++) { 
+      
+      Int_t left   = fADCT[iCol-1][iT]; 
+      Int_t center = fADCT[iCol][iT];
+      Int_t right  = fADCT[iCol+1][iT];  
+
+      Int_t sum = left + center + right;            // cluster charge sum
+      qsumAlu.AssignFormatted(sum);    
+      qsumAlu = qsumAlu;                        // size-checking; redundant
+      qsum[iCol] = qsumAlu.GetValue(); 
+      
+      //hit detection and masking
+      if(center>=left){
+       if(center>right){
+         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
+           mark |= 1;                              // marker
+         }
+       }
+      }
+      mark = mark<<1;                
+    }
+    mark = mark>>1;
+
+       
+    // get selection of 6 adc's and sort,starting with greatest values
+
+    //read three from right side and sort (primitive sorting algorithm)
+    Int_t i = 0; // adc number
+    Int_t j = 1; // selection number
+    while(i<fNADC-2 && j<=3){
+      i = i + 1;
+      if((mark>>(i-1)) & 1 == 1) {
+       selection[j].iadc  = fNADC-1-i;
+       selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
+       
+       // insert into sorted list
+       listLeft = &selection[0];
+       list = listLeft->next;
+       
+       if(list!=NULL) {
+         while((list->next != NULL) && (selection[j].value <= list->value)){
+           listLeft = list;
+           list = list->next;
+         }
+         
+         if(selection[j].value<=list->value){
+           selection[j].next = list->next;
+           list->next = &selection[j];
+         }
+         else {
+           listLeft->next = &selection[j];
+           selection[j].next = list;
+         }
+       }
+       else{
+         listLeft->next = &selection[j];
+         selection[j].next = list;
+       }
+       
+       j = j + 1;
+      }
+    }
+
+
+    // read three from left side
+    Int_t k = fNADC-2;
+    while(k>i && j<=6) {
+      if((mark>>(k-1)) & 1 == 1) {
+       selection[j].iadc  = fNADC-1-k;
+       selection[j].value = qsum[fNADC-1-k]>>6;
+       
+       listLeft = &selection[0];
+       list = listLeft->next;
+       
+       if(list!=NULL){
+         while((list->next != NULL) && (selection[j].value <= list->value)){
+           listLeft = list;
+           list = list->next;
+         }
+       
+         if(selection[j].value<=list->value){
+           selection[j].next = list->next;
+           list->next = &selection[j];
+         }
+         else {
+           listLeft->next = &selection[j];
+           selection[j].next = list;
+         }
+       }
+       else{
+         listLeft->next = &selection[j];
+         selection[j].next = list;
+       }
+
+       j = j + 1;
+      }
+      k = k - 1;
+    }
+
+    // get the four with greatest charge-sum
+    list = &selection[0];
+    for(i = 0; i<4; i++){
+      if(list->next == NULL) continue;
+      list = list->next;
+      if(list->iadc == -1) continue;
+      Int_t adc = list->iadc;                              // channel number with selected hit
+      
+      // the following arrays contain the four chosen channels in 1 time-bin
+      portChannel[i]   = adc; 
+      clusterCharge[i] = qsum[adc];
+      leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
+      centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];           
+      rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];         
+    }
+
+    // arithmetic unit
+    
+    // cluster verification
+    if(!bVBY){
+      for(i = 0; i<4; i++){
+       Int_t lr = leftCharge[i]*rightCharge[i]*1024;
+       Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
+       if (lr>=cc){
+         portChannel[i]   = -1;                                 // set to invalid address 
+         clusterCharge[i] = 0;
+       }
+      }
+    }
+
+    // fit-sums of valid channels
+    // local hit position
+    for(i = 0; i<4; i++){
+      if (centerCharge[i] ==  0) {
+       portChannel[i] = -1; 
+      }// prevent division by 0
+      
+      if (portChannel[i]  == -1) continue;
+      
+      Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
+       
+      Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
+      dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
+      dCOGAlu.AssignDouble(dCOG);
+      Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
+            
+      dCOG = dCOG/2;
+      yrawAlu.AssignDouble(dCOG);
+      Int_t iCOG = yrawAlu.GetValue();
+      Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
+      yrawAlu.AssignFormatted(y);               // 0<y<1           
+      yAlu  = yrawAlu;                        // convert to 16 past-comma bits
+      
+      if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
+      xAlu.AssignInt(iT);                       // buffer width of 5 bits 
+      
+
+      xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x       
+      
+      yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
+   
+      xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
+                  
+      Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
+
+      // calculate fit-sums recursively
+      // interpretation of their bit-length is given as comment
+      
+      // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
+
+      XAlu.AssignFormatted(X[adc]);
+      XAlu = XAlu + xAlu;                   // buffer width of 9 bits 
+      X[adc] = XAlu.GetValue();
+             
+      XXAlu.AssignFormatted(XX[adc]);
+      XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits    
+      XX[adc] = XXAlu.GetValue();
+
+      if (Y[adc] < 0) {
+       YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
+       YAlu.SetSign(-1);
+      }
+      else {
+       YAlu.AssignFormatted(Y[adc]);
+       YAlu.SetSign(1);
+      }
+       
+      YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);     
+      Y[adc] = YAlu.GetSignedValue();
+            
+      YYAlu.AssignFormatted(YY[adc]);
+      YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma) 
+      YY[adc] = YYAlu.GetValue();
+           
+      if (XY[adc] < 0) {
+       XYAlu.AssignFormatted(-XY[adc]);
+       XYAlu.SetSign(-1);
+      }
+      else {
+       XYAlu.AssignFormatted(XY[adc]);
+       XYAlu.SetSign(1);
+      }
+
+      XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma) 
+      XY[adc] = XYAlu.GetSignedValue();
+            
+      N[adc]  = N[adc] + 1;
+   
+
+      // accumulated charge
+      qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
+      qtruncAlu = qsumAlu;
+
+      if(iT>=tQS0 && iT<=tQE0){
+       QT0Alu.AssignFormatted(QT0[adc]);
+       QT0Alu = QT0Alu + qtruncAlu;
+       QT0[adc] = QT0Alu.GetValue();
+       //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
+      }
+      
+      if(iT>=tQS1 && iT<=tQE1){
+       QT1Alu.AssignFormatted(QT1[adc]);
+       QT1Alu = QT1Alu + qtruncAlu;
+       QT1[adc] = QT1Alu.GetValue();
+       //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
+      }
+    }// i
+      
+    // remapping is done!!
+     
+  }//iT
+  
+    
+  // tracklet-assembly
+  
+  // put into AliTRDfeeParam and take care that values are in proper range
+  const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??) 
+  const Int_t cTCT = 8;      // joint number of hits;     8<=TCT<=31
+  
+  Int_t mPair   = 0;         // marker for possible tracklet pairs
+  Int_t* hitSum = new Int_t[fNADC-3];
+  // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18; 
+  
+  // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
+  for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
+    hitSum[iCol] = N[iCol] + N[iCol+1];
+    if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
+       mPair |= 1;         // mark as possible channel-pair
+     
+    }
+    mPair = mPair<<1;
+  }
+  mPair = mPair>>1;
+  
+  List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
+                                                 // selectPair[18] is starting list-element just for pointing
+  for(Int_t k = 0; k<fNADC-2; k++){
+      selectPair[k].next       = NULL;
+      selectPair[k].iadc       =   -1;           // invalid adc
+      selectPair[k].value      =    0;
+   
+    }
+
+ list = NULL;
+ listLeft = NULL;
+  
+  // read marker and sort according to hit-sum
+  
+  Int_t adcL  = 0;            // left adc-channel-number (remapped)
+  Int_t selNr = 0;            // current number in list
+  
+  // insert marked channels into list and sort according to hit-sum
+  while(adcL < fNADC-3 && selNr < fNADC-3){
+     
+    if((mPair>>((fNADC-4)-(adcL))) & 1 == 1) {
+      selectPair[selNr].iadc  = adcL;
+      selectPair[selNr].value = hitSum[adcL];   
+      
+      listLeft = &selectPair[fNADC-3];
+      list = listLeft->next;
+       
+      if(list!=NULL) {
+       while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
+         listLeft = list;
+         list = list->next;
+       }
+       
+       if(selectPair[selNr].value <= list->value){
+         selectPair[selNr].next = list->next;
+         list->next = &selectPair[selNr];
+       }
+       else {
+         listLeft->next = &selectPair[selNr];
+         selectPair[selNr].next = list;
+       }
+       
+      }
+      else{
+       listLeft->next = &selectPair[selNr];
+       selectPair[selNr].next = list;
+      }
+      
+      selNr = selNr + 1;
+    }
+    adcL = adcL + 1;
+  }
+  
+  //select up to 4 channels with maximum number of hits
+  Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
+  Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
+  list = &selectPair[fNADC-3];
+  
+  for (Int_t i = 0; i<4; i++) {
+    if(list->next == NULL) continue;
+    list = list->next;
+    if(list->iadc == -1) continue;
+    lpairChannel[i] = list->iadc;        // channel number with selected hit
+    rpairChannel[i] = lpairChannel[i]+1;
+  }
+  
+  // avoid submission of double tracklets  
+  for (Int_t i = 3; i>0; i--) {
+    for (Int_t j = i-1; j>-1; j--) {
+      if(lpairChannel[i] == rpairChannel[j]) {
+       lpairChannel[i] = -1;
+       rpairChannel[i] = -1;
+       break;
+      }
+      if(rpairChannel[i] == lpairChannel[j]) {
+       lpairChannel[i] = -1;
+       rpairChannel[i] = -1;
+       break;
+      }
+    }
+  }
+  
+  // merging of the fit-sums of the remainig channels
+  // assume same data-word-width as for fit-sums for 1 channel
+  // relative scales!
+  Int_t mADC[4];                      
+  Int_t mN[4];
+  Int_t mQT0[4];
+  Int_t mQT1[4];
+  Int_t mX[4];
+  Int_t mXX[4];
+  Int_t mY[4];
+  Int_t mYY[4];
+  Int_t mXY[4];
+  Int_t mOffset[4];
+  Int_t mSlope[4];
+  Int_t mMeanCharge[4]; 
+  Int_t inverseN = 0;
+  Double_t invN = 0;
+  for (Int_t i = 0; i<4; i++){
+    mADC[i] = -1;                        // set to invalid number
+    mN[i]   =  0;
+    mQT0[i] =  0;
+    mQT1[i] =  0;
+    mX[i]   =  0;
+    mXX[i]  =  0;
+    mY[i]   =  0;
+    mYY[i]  =  0;
+    mXY[i]  =  0;
+    mOffset[i] = 0;
+    mSlope[i]  = 0;
+    mMeanCharge[i] = 0;
+  }
+  
+
+  YAlu.AssignInt(1);
+  Int_t wpad  = YAlu.GetValue();       // 1 with 8 past-comma bits
+  for (Int_t i = 0; i<4; i++){
+    
+    mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
+                                         // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
+    Int_t madc = mADC[i];
+    if (madc == -1) continue;
+    mN[i]    = hitSum[madc];
+  
+    // 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
+    if (N[madc+1] == 0) {
+       mQT0[i] = QT0[madc];
+       mQT1[i] = QT1[madc];
+       
+    }
+    else {
+
+       // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
+       
+       mQT0[i]   = QT0[madc] + QT0[madc+1];
+       QT0Alu.AssignFormatted(mQT0[i]);   
+       QT0Alu  = QT0Alu;                // size-check
+       mQT0[i]   = QT0Alu.GetValue();     // write back
+       
+       mQT1[i]   = QT1[madc] + QT1[madc+1];
+       QT1Alu.AssignFormatted(mQT1[i]);
+       QT1Alu  = QT1Alu;
+       mQT1[i]   = QT1Alu.GetValue();
+    }
+    
+    // calculate the mean charge in adc values; later to be replaced by electron likelihood
+    mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
+    mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
+    // 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
+    invN = 1.0/(mN[i]);
+    inverseNAlu.AssignDouble(invN);
+    inverseN = inverseNAlu.GetValue();
+    mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
+    TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
+    TotalChargeAlu = TotalChargeAlu;
+    MeanChargeAlu = TotalChargeAlu;
+    mMeanCharge[i] = MeanChargeAlu.GetValue();
+    
+    if (N[madc+1] == 0) {
+       mX[i]     =   X[madc];
+       mXX[i]    =  XX[madc];
+       mY[i]     =   Y[madc];
+       mXY[i]    =  XY[madc];
+       mYY[i]    =  YY[madc];
+    }
+    else {
+       
+       mX[i]     =   X[madc] +  X[madc+1];
+       XAlu.AssignFormatted(mX[i]);
+       XAlu    = XAlu;
+       mX[i]     = XAlu.GetValue();
+       
+       mXX[i]    =  XX[madc] + XX[madc+1];
+       XXAlu.AssignFormatted(mXX[i]);
+       XXAlu   = XXAlu;
+       mXX[i]    = XXAlu.GetValue();
+    
+       mY[i]     =   Y[madc] + Y[madc+1] + wpad;
+       if (mY[i] < 0) {
+           YAlu.AssignFormatted(-mY[i]);
+           YAlu.SetSign(-1);
+       }
+       else {
+           YAlu.AssignFormatted(mY[i]);
+           YAlu.SetSign(1);
+       }
+       YAlu    = YAlu;
+       mY[i]     = YAlu.GetSignedValue();
+       
+       mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*wpad;
+       
+       if (mXY[i] < 0) {
+           XYAlu.AssignFormatted(-mXY[i]);
+           XYAlu.SetSign(-1);
+       }
+       else {
+           XYAlu.AssignFormatted(mXY[i]);
+           XYAlu.SetSign(1);
+       }
+       XYAlu   = XYAlu;
+       mXY[i]    = XYAlu.GetSignedValue();
+    
+       mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*wpad + wpad*wpad;
+       if (mYY[i] < 0) {
+           YYAlu.AssignFormatted(-mYY[i]);
+           YYAlu.SetSign(-1);
+       }
+       else {
+           YYAlu.AssignFormatted(mYY[i]);
+           YYAlu.SetSign(1);
+       }
+       
+       YYAlu   = YYAlu;
+       mYY[i]    = YYAlu.GetSignedValue();
+    }
+  
+  }
+    
+  // calculation of offset and slope from the merged fit-sums; 
+  // YY is needed for some error measure only; still to be done
+  // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
+  
+  // which formats should be chosen?
+  AliTRDtrapAlu denomAlu;
+  denomAlu.Init(20,8);       
+  AliTRDtrapAlu numAlu;
+  numAlu.Init(20,8);     
+  // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
+  // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
+
+  for (Int_t i = 0; i<4; i++) {
+    if (mADC[i] == -1) continue;
+      
+    Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
+    if (num0 < 0) {
+      denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
+      denomAlu.SetSign(-1);
+    }
+    else {
+      denomAlu.AssignInt(num0);
+      denomAlu.SetSign(1);
+    }
+    
+    Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
+    if (num1 < 0) {
+      numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
+      numAlu.SetSign(-1);
+    }
+    else {
+      numAlu.AssignFormatted(num1);
+      numAlu.SetSign(1);
+    }
+    numAlu    = numAlu/denomAlu;
+    mSlope[i]   = numAlu.GetSignedValue();
+   
+    Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
+   
+    if (num2 < 0) {
+      numAlu.AssignFormatted(-num2);
+      numAlu.SetSign(-1);
+    }
+    else {
+      numAlu.AssignFormatted(num2);
+      numAlu.SetSign(1);
+    }
+   
+    numAlu    = numAlu/denomAlu;
+   
+    
+    mOffset[i]  = numAlu.GetSignedValue();
+    numAlu.SetSign(1);
+    denomAlu.SetSign(1);
+       
+                                 
+    //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)             
+    numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
+    mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values 
+    
+    // up to here: adc-mapping according to TRAP and in line with pad-col mapping
+    // reverse adc-counting to be again in line with the online mapping
+    mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
+    mADC[i]     = mADC[i] + 2; 
+    // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
+  }
+
+  // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected; 
+  // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
+  
+  // transform parameters to the local coordinate-system of a stack (used by GTU)
+  AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
+  
+  Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
+  //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
+  
+  // difference between width of inner and outer pads of a row is not accounted for;
+  
+  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 
+  Double_t eCharge  = 0.3;                            // unit charge in (GeV/c)/m*T
+  Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
+  
+  Double_t granularityOffset = 0.160;                // granularity for offset in mm
+  Double_t granularitySlope  = 0.140;                // granularity for slope  in mm     
+    
+  // get the coordinates in SM-system; parameters: 
+  
+  Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5 
+  // zPos is position of pad-borders;
+  Double_t zOffset = 0.0;
+  if ( fRow == 0 || fRow == 15 ) {
+      zOffset = padPlane->GetLengthOPad();
+  }
+  else {
+      zOffset = padPlane->GetLengthIPad();
+  }
+  zOffset = (-1.0) * zOffset/2.0;
+  // turn zPos to be z-coordinate at middle of pad-row
+  zPos = zPos + zOffset;
+
+      
+  Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
+  Int_t    icol        =  0;                                 // column-number of adc-channel
+  Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
+  Double_t dx          = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
+  Double_t vdrift      = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector
+  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
+  //Double_t lorAngle   =  7.0;                               // Lorentz-angle in degrees
+  Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
+  Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
+  //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
+
+  Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
+  Double_t slopeMin[4];                            // local limits for the deflection
+  Double_t slopeMax[4];
+  Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
+  Int_t   mslopeMax[4];
+
+
+  //x coord. of upper side of drift chambers in local SM-system (in mm)
+  //obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
+  //the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
+  switch(fLayer) 
+    {
+    case 0: 
+      xPos = 3003.0;
+      break;
+    case 1:
+      xPos = 3129.0;
+      break;
+    case 2:
+      xPos = 3255.0;
+      break;
+    case 3:
+      xPos = 3381.0;
+      break;
+    case 4:
+      xPos = 3507.0;
+      break;
+    case 5:
+      xPos = 3633.0;
+      break;
+    }
+  // calculation of offset-correction n: 
+
+  Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));  
+  nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
+  if (nCorrectOffset < 0) {
+    numAlu.AssignInt(-nCorrectOffset);
+    numAlu.SetSign(-1);
+  }
+  else {
+    numAlu.AssignInt(nCorrectOffset);
+    numAlu.SetSign(1);
+  }
+  nCorrectOffset = numAlu.GetSignedValue();         
+  Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
+  // calculation of slope-correction
+
+  // this is only true for tracks coming (approx.) from primary vertex
+  Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
+  // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
+  // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
+  
+  if (cCorrectSlope < 0) {
+      numAlu.AssignDouble(-cCorrectSlope);
+      numAlu.SetSign(-1);
+  }
+  else {
+      numAlu.AssignDouble(cCorrectSlope);
+      numAlu.SetSign(1);
+  }
+  cCorrectSlope = numAlu.GetSignedValue();
+  // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
+  // different pad-width of outer pads of a pad-plane not taken into account
+  // tFS: upper plane of drift-volume (not amplification region); this choice is important for offset
+  // tFE: does !!not!! need to correspond to lower plane of drift-volume; TR hits can be cut;
+  
+  Double_t mCorrectSlope = ((Double_t)(fNTimeBin-tFS))*padWidthI/granularitySlope;  // >= 0.0
+
+  AliTRDtrapAlu correctAlu;
+  correctAlu.Init(20,8);
+  
+  AliTRDtrapAlu offsetAlu;
+  offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
+  
+  AliTRDtrapAlu slopeAlu;
+  slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
+
+  for (Int_t i = 0; i<4; i++) {
+    
+    if (mADC[i] == -1) continue;
+    
+    icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
+    yPos[i]   = (padPlane->GetColPos(icol))*10.0;
+    
+    
+    // offset:
+    
+    correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
+    mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
+    mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset)); 
+    mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
+    
+    if (mOffset[i] < 0) {
+      numAlu.AssignFormatted(-mOffset[i]);
+      numAlu.SetSign(-1);
+    }
+    else {
+      numAlu.AssignFormatted(mOffset[i]);
+      numAlu.SetSign(1);
+    }
+
+    offsetAlu = numAlu; 
+    mOffset[i]  = offsetAlu.GetSignedValue();  
+
+    
+    // slope:
+    
+    correctAlu.AssignDouble(mCorrectSlope);
+    mCorrectSlope = correctAlu.GetValueWhole();
+    
+    mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
+
+    if (mSlope[i] < 0) {
+      numAlu.AssignFormatted(-mSlope[i]);
+      numAlu.SetSign(-1);
+    }
+    else {
+      numAlu.AssignFormatted(mSlope[i]);
+      numAlu.SetSign(1);
+    }
+
+    slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
+    mSlope[i]   = slopeAlu.GetSignedValue(); 
+       
+    // local (LTU) limits for the deflection 
+    // ATan returns angles in radian
+    alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
+    slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
+    slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
+    
+    if (slopeMin[i] < 0) {
+      slopeAlu.AssignDouble(-slopeMin[i]);
+      slopeAlu.SetSign(-1);
+    }
+    else { 
+      slopeAlu.AssignDouble(slopeMin[i]);
+      slopeAlu.SetSign(1);
+    }
+    mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
+   
+    if (slopeMax[i] < 0) {
+      slopeAlu.AssignDouble(-slopeMax[i]);
+      slopeAlu.SetSign(-1);
+    }
+    else {
+      slopeAlu.AssignDouble(slopeMax[i]);
+      slopeAlu.SetSign(1);
+    }
+    mslopeMax[i] = slopeAlu.GetSignedValue();
+  }
+
+  // suppress submission of tracks with low stiffness
+  // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber) 
+
+  // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
+  // up to now they are sorted according to maximum hit sum
+  // is the sorting really done in the TRAP-chip?
+  
+  Int_t order[4] = {-1,-1,-1,-1};
+  Int_t wordnr = 0;   // number of tracklet-words
+  
+  for(Int_t j = 0; j < fMaxTracklets; j++) {
+      //if( mADC[j] == -1) continue; 
+      if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
+      wordnr++;
+      if( wordnr-1 == 0) {
+         order[0] = j;
+         continue;
+      }
+      // wordnr-1>0, wordnr-1<4
+      order[wordnr-1] = j;
+      for( Int_t k = 0; k < wordnr-1; k++) {
+         if( mOffset[j] < mOffset[order[k]] ) {
+             for( Int_t l = wordnr-1; l > k; l-- ) {
+                 order[l] = order[l-1];
+             }
+             order[k] = j;
+             break;
+         }
+         
+      }
+  }
+        
+  // fill the bit-words in ascending order and without gaps
+  UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
+  for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
+      //Bool_t rem1 = kTRUE;
+    
+    Int_t i = order[j];
+    bitWord[j] =   0; // invalid bit-word (bit-word is 2-complement and therefore without sign)
+    //if( mADC[i] == -1) continue; 
+    ////if( (mADC[i] == -1) || (mSlope[i] < mslopeMin[i]) || (mSlope[i] > mslopeMax[i])) continue; //don't transmit bit word 
+    bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
+    //printf("\n");
+    
+    /*
+    // pad position
+    if(mOffset[i] < 0) {
+      rem1 = kFALSE;   // don't remove the first 1
+      //printf("1");
+      for(Int_t iBit = 1; iBit < 13; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (1-((-mOffset[i])>>(12-iBit))&1);
+       //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
+      }
+    }
+    else {
+      bitWord[j]   |= 0; 
+      //printf("0");
+      for(Int_t iBit = 1; iBit < 13; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (mOffset[i]>>(12-iBit))&1;
+       //printf("%d",(mOffset[i]>>(12-iBit))&1);
+      }
+    }
+        
+    // deflection length
+    bitWord[j] = bitWord[j]<<1;
+    if(mSlope[i] < 0) {
+      bitWord[j]   |= 1;
+      //printf("1");
+      for(Int_t iBit = 1; iBit < 7; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (1-((-mSlope[i])>>(6-iBit))&1);
+       //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
+      }
+    }
+    else {
+      bitWord[j]   |= 0;
+      //printf("0");
+      for(Int_t iBit = 1; iBit < 7; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (mSlope[i]>>(6-iBit))&1;
+       //printf("%d",(mSlope[i]>>(6-iBit))&1);
+      }
+    }
+
+    // pad row
+    for(Int_t iBit = 0; iBit < 4; iBit++) {
+      bitWord[j]  = bitWord[j]<<1;
+      bitWord[j] |= (fRow>>(3-iBit))&1;
+      //printf("%d", (fRow>>(3-iBit))&1);
+    }
+    
+    // electron probability (currently not implemented; the mean charge is just scaled)
+    for(Int_t iBit = 0; iBit < 8; iBit++) {
+      bitWord[j]  = bitWord[j]<<1;
+      bitWord[j] |= (mMeanCharge[i]>>(7-iBit))&1;               
+      //printf("0");
+    }
+    
+
+    if (rem1 == kTRUE) {
+      bitWord[j] = bitWord[j] - (1<<31);
+      }*/
+
+
+    /*printf("mean charge: %d\n",mMeanCharge[i]);
+    printf("row: %d\n",fRow);
+    printf("slope: %d\n",mSlope[i]);
+    printf("pad position: %d\n",mOffset[i]);
+    printf("channel: %d\n",mADC[i]);*/
+
+    // electron probability (currently not implemented; the mean charge is just scaled)
+    for(Int_t iBit = 0; iBit < 8; iBit++) {
+      bitWord[j]  = bitWord[j]<<1;
+      bitWord[j] |= (mMeanCharge[i]>>(7-iBit))&1;               
+      //printf("0");
+    }
+    
+    // pad row
+    for(Int_t iBit = 0; iBit < 4; iBit++) {
+      bitWord[j]  = bitWord[j]<<1;
+      bitWord[j] |= (fRow>>(3-iBit))&1;
+      //printf("%d", (fRow>>(3-iBit))&1);
+    }
+    
+    // deflection length
+    bitWord[j] = bitWord[j]<<1;
+    if(mSlope[i] < 0) {
+      bitWord[j]   |= 1;
+      //printf("1");
+      for(Int_t iBit = 1; iBit < 7; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (1-((-mSlope[i])>>(6-iBit))&1);
+       //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
+      }
+    }
+    else {
+      bitWord[j]   |= 0;
+      //printf("0");
+      for(Int_t iBit = 1; iBit < 7; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (mSlope[i]>>(6-iBit))&1;
+       //printf("%d",(mSlope[i]>>(6-iBit))&1);
+      }
+    }
+
+    // pad position
+    bitWord[j] = bitWord[j]<<1;
+    if(mOffset[i] < 0) {
+       bitWord[j]   |= 1;
+       //printf("1");
+       for(Int_t iBit = 1; iBit < 13; iBit++) {
+           bitWord[j]  = bitWord[j]<<1;
+           bitWord[j] |= (1-((-mOffset[i])>>(12-iBit))&1);
+           //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
+       }
+    }
+    else {
+      bitWord[j]   |= 0; 
+      //printf("0");
+      for(Int_t iBit = 1; iBit < 13; iBit++) {
+       bitWord[j]  = bitWord[j]<<1;
+       bitWord[j] |= (mOffset[i]>>(12-iBit))&1;
+       //printf("%d",(mOffset[i]>>(12-iBit))&1);
+      }
+    }
+
+
+   
+    
+    //printf("bitWord: %u\n",bitWord[j]);
+    //printf("adc: %d\n",mADC[i]);
+    fMCMT[j] = bitWord[j];
+  }
+    
+  //printf("\n");
+
+  
+  delete [] qsum;
+  delete [] ieffped;
+
+  delete [] X;
+  delete [] XX;
+  delete [] Y;
+  delete [] YY;
+  delete [] XY;
+  delete [] N;
+  delete [] QT0;
+  delete [] QT1;
+
+  delete [] hitSum;
+  delete [] selectPair;
+
+  delete padPlane;
+
+
+
+/*
+
+  // output-part; creates some dump trees; output should not be organized inside the AliTRDmcmSim-class
+  
+  // structure: in system directory "./TRD_Tracklet" a root-file called "TRD_readout_tree.root" is stored with subdirectories SMxx/sx;
+  // in each of these subdirectories 6 trees according to layers are saved, called lx; 
+  // whenever a mcm of that layer had a bit-word!=0, a branch containing an array with 4 (possibly some valued 0) elements is added;
+  // branch-name: mcmxxxwd; 
+  // another branch contains the channel-number (mcmxxxch)
+  
+  
+  AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
+  AliLog::SetFileOutput("../log/tracklet.log");
+  
+  UInt_t* trackletWord;
+  Int_t*  adcChannel;
+
+  Int_t u = 0;
+    
+  // testing for wordnr in order to speed up the simulation
+  
+  if (wordnr == 0 && fNextEvent == 0) {
+    return;
+  }
+
+   
+  Int_t mcmNr = fRobPos * (fGeo->MCMmax()) + fMcmPos;
+  
+  Char_t* SMName    = new Char_t[4];
+  Char_t* stackName = new Char_t[2];
+  Char_t* layerName = new Char_t[2];
+
+  Char_t* treeName  = new Char_t[2];
+  Char_t* treeTitle = new Char_t[8];
+  Char_t* branchNameWd = new Char_t[8]; // mcmxxxwd bit word
+  Char_t* branchNameCh = new Char_t[8]; // mcmxxxch channel
+   
+  Char_t* dirName = NULL;
+  Char_t* treeFile  = NULL; 
+  Char_t* evFile = NULL;
+  Char_t* curDir = new Char_t[255];
+  
+  for (Int_t i = 0; i<255; i++) {
+    curDir[i]='n';
+  }
+  sprintf(curDir,"%s",gSystem->BaseName(gSystem->WorkingDirectory()));
+  Int_t rawEvent = 0;
+  Int_t nrPos = 3;
+  
+
+  while(curDir[nrPos]!='n'){
+    
+   
+    switch(curDir[nrPos]) {
+    case '0':
+      rawEvent = rawEvent*10 + 0;
+      break;
+    case '1':
+      rawEvent = rawEvent*10 + 1;
+      break;
+    case '2':
+      rawEvent = rawEvent*10 + 2;
+      break;
+    case '3':
+      rawEvent = rawEvent*10 + 3;
+      break;
+    case '4':
+      rawEvent = rawEvent*10 + 4;
+      break;
+    case '5':
+      rawEvent = rawEvent*10 + 5;
+      break;
+    case '6':
+      rawEvent = rawEvent*10 + 6;
+      break;
+    case '7':
+      rawEvent = rawEvent*10 + 7;
+      break;
+    case '8':
+      rawEvent = rawEvent*10 + 8;
+      break;
+    case '9':
+      rawEvent = rawEvent*10 + 9;
+      break;
+   
+    }
+    nrPos = nrPos + 1; 
+  }
+  delete [] curDir;
+    
+  if (!gSystem->ChangeDirectory("../TRD_Tracklet")) {
+    gSystem->MakeDirectory("../TRD_Tracklet");
+    gSystem->ChangeDirectory("../TRD_Tracklet");
+  }
+  
+  TFile *f = new TFile("TRD_readout_tree.root","update");
+  TTree *tree          = NULL;
+  TBranch *branch      = NULL;
+  TBranch *branchCh   = NULL; 
+   
+  Int_t iEventNr = 0; 
+  Int_t dignr = 10; // nr of digits of a integer
+  Int_t space = 1;  // additional char-space
+  
+  
+  evFile = new Char_t[2+space];
+  sprintf(evFile,"ev%d",iEventNr);
+
+  
+  while(f->cd(evFile)){
+    iEventNr = iEventNr + 1;
+    if (iEventNr/dignr > 0) {
+      dignr = dignr * 10;
+      space = space + 1;
+    }
+    delete [] evFile;
+    evFile = NULL;
+    evFile = new Char_t[2+space];
+    sprintf(evFile,"ev%d",iEventNr); 
+  }
+  
+  if(iEventNr == rawEvent) { fNextEvent = 1; } // new event
+   
+  if (fNextEvent == 1) {
+    fNextEvent = 0;
+    // turn to head directory
+    f->mkdir(evFile);
+    f->cd(evFile);
+    // create all subdirectories and trees in case of new event
+   
+     
+    for (Int_t iSector = 0; iSector < 18; iSector++) {
+  
+      if (iSector < 10) {
+       sprintf(SMName,"SM0%d",iSector);
+      }
+      else {
+       sprintf(SMName,"SM%d",iSector);
+      }
+
+      for (Int_t iStack = 0; iStack < 5; iStack++) {
+       sprintf(stackName,"s%d",iStack);
+       
+       f->cd(evFile);
+       if (iStack == 0) {
+         gDirectory->mkdir(SMName);
+       }
+       gDirectory->cd(SMName);
+       gDirectory->mkdir(stackName);
+       gDirectory->cd(stackName);
+       
+       for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+         sprintf(layerName,"l%d",iLayer);
+         sprintf(treeName,"%s",layerName);
+         sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
+         tree = new TTree(treeName,treeTitle);
+         tree->Write("",TObject::kOverwrite);
+         delete tree;
+         tree = NULL;
+       }
+      }
+    }
+   
+    
+  }
+
+    
+  else {
+    iEventNr = iEventNr - 1;
+    dignr = dignr/10;
+    if (iEventNr/dignr == 0) space = space - 1;
+    delete [] evFile;
+    evFile = NULL;
+    evFile = new Char_t[2+space];
+    sprintf(evFile,"ev%d",iEventNr);    
+  }
+  
+  if (wordnr == 0) {
+    delete [] SMName;
+    delete [] stackName;
+    delete [] layerName;
+    delete [] treeName;
+    delete [] treeTitle;
+    delete [] branchNameWd;
+    delete [] branchNameCh;
+    delete [] evFile;
+    f->Close();
+    dirName   = new Char_t[6+space];
+    sprintf(dirName,"../raw%d",iEventNr);
+    gSystem->ChangeDirectory(dirName);
+    delete [] dirName;
+    return;
+  }
+  
+  dirName   = new Char_t[6+space];
+  sprintf(dirName,"../raw%d",iEventNr);
+  f->cd(evFile);
+  if (fSector < 10) {
+    sprintf(SMName,"SM0%d",fSector);
+  }
+  else {
+    sprintf(SMName,"SM%d",fSector);
+  }
+  sprintf(stackName,"s%d",fStack);
+  sprintf(layerName,"l%d",fLayer);
+  sprintf(treeName,"%s",layerName);
+  sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
+  treeFile = new Char_t[13+space];
+  sprintf(treeFile,"%s/%s/%s/%s",evFile,SMName,stackName,treeName);
+  delete [] evFile;
+  evFile = NULL;
+  gDirectory->cd(SMName);
+  gDirectory->cd(stackName);
+  tree = (TTree*)f->Get(treeFile);
+  delete [] treeFile;
+  treeFile = NULL;
+
+
+  //make branch with number of words and fill
+
+  if (mcmNr < 10) {
+    sprintf(branchNameWd,"mcm00%dwd",mcmNr);
+    sprintf(branchNameCh,"mcm00%dch",mcmNr);
+  }
+  else {
+    if (mcmNr < 100) {
+      sprintf(branchNameWd,"mcm0%dwd",mcmNr); 
+      sprintf(branchNameCh,"mcm0%dch",mcmNr);
+    }
+    else {
+      sprintf(branchNameWd,"mcm%dwd",mcmNr); 
+      sprintf(branchNameCh,"mcm%dch",mcmNr);
+    }
+  }
+
+      
+  
+  // fill the tracklet word; here wordnr > 0
+
+  trackletWord = new UInt_t[fMaxTracklets];
+  adcChannel   = new Int_t[fMaxTracklets];
+
+  for (Int_t j = 0; j < fMaxTracklets; j++) {
+      Int_t i = order[j];
+      trackletWord[j] = 0;
+      adcChannel[j] = -1;
+      if (bitWord[j]!=0) {
+         trackletWord[u] = bitWord[j];
+         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
+         
+         //fMCMT[u] = bitWord[j];
+         u = u + 1;
+      }
+  }
+  
+  branch = tree->GetBranch(branchNameWd);
+  if(!branch) {
+    //make branch and fill
+    branch = tree->Branch(branchNameWd,trackletWord,"trackletWord[4]/i"); // 32 bit unsigned integer
+    branch->Fill();
+  }
+
+  branchCh = tree->GetBranch(branchNameCh);
+  if(!branchCh) {
+    //make branch and fill
+    branchCh = tree->Branch(branchNameCh,adcChannel,"adcChannel[4]/i"); // 32 bit unsigned integer
+    branchCh->Fill();
+  }
+
+  tree->Write("",TObject::kOverwrite);
+
+  delete [] SMName;
+  delete [] stackName;
+  delete [] layerName;
+  delete [] treeName;
+  delete [] treeTitle;
+  delete [] branchNameWd;
+  delete [] branchNameCh;
+  delete [] trackletWord;
+  delete [] adcChannel;
+  
+  f->Close();
+  gSystem->ChangeDirectory(dirName);
+  delete [] dirName;
+  
+*/
+
+  // to be done:
+  // error measure for quality of fit (not necessarily needed for the trigger)
+  // cluster quality threshold (not yet set)
+  // electron probability
+  
+
+   
+}
+
+
+
+
+
+
index e198d86..de02901 100644 (file)
@@ -17,6 +17,8 @@ class AliTRDfeeParam;
 class AliTRDSimParam;
 class AliTRDcalibDB;
 class AliTRDgeometry;
+class AliTRDtrapAlu;
+class AliTRDpadPlane;
 
 class AliTRDmcmSim : public TObject {
 
@@ -29,6 +31,7 @@ class AliTRDmcmSim : public TObject {
 
   virtual void      Copy(TObject &m) const;
 
+  // void      Init( Int_t cha, Int_t rob, Int_t mcm, Bool_t newEvent );   // Initialize MCM by the position parameters
           void      Init( Int_t cha, Int_t rob, Int_t mcm );   // Initialize MCM by the position parameters
           void      SetData(Int_t iadc, Int_t *adc);           // Set ADC data with array 
           void      SetData(Int_t iadc, Int_t it, Int_t adc ); // Set ADC data
@@ -41,14 +44,21 @@ class AliTRDmcmSim : public TObject {
          Int_t     GetCol( Int_t iadc );                      // Get corresponding column (0-143) from for ADC channel iadc = [0:20]
          // for the ADC/Col mapping, see: http://wiki.kip.uni-heidelberg.de/ti/TRD/index.php/Image:ROB_MCM_numbering.pdf
 
+         Int_t*    GetPosLUT();
+
          Int_t     ProduceRawStream( UInt_t *buf, Int_t bufsize ); // Produce raw data stream from this MCM
+         Int_t     ProduceTrackletStream( UInt_t *buf, Int_t bufsize ); // produce the tracklet stream for this MCM
          void      Filter();                                  // Apply digital filters for existing data
          void      ZSMapping();                               // Do ZS mapping for existing data
          void      DumpData( char *f, char *target );         // Dump data stored (only for debugging)
+         void      Tracklet();
+         void      SetPosLUT();
 
  protected:
 
           Bool_t    fInitialized;                       // Status whether the class is initialized or not
+         Bool_t    fNextEvent;                         // 0, if new event; 1 if not
+         Int_t     fMaxTracklets;                      // maximum number of tracklet-words submitted per mcm **new**
          Int_t     fChaId;                             // Chamber ID
          Int_t     fSector;                            // Sector number of the supermodule
          Int_t     fStack;                             // Chamber stack ID
@@ -58,11 +68,18 @@ class AliTRDmcmSim : public TObject {
          Int_t     fNADC;                              // Number of ADC (usually 21)
          Int_t     fNTimeBin;                          // Number of Timebins (variable)
           Int_t     fRow;                               // Pad row number (0-11 or 0-15) of the MCM on chamber
-          Int_t   **fADCR;                              // Array with MCM ADC values (Raw)
+         Int_t   **fADCR;                              // Array with MCM ADC values (Raw)
           Int_t   **fADCF;                              // Array with MCM ADC values (Filtered)
-          Int_t   **fZSM;                               // Zero suppression map
+         Int_t   **fADCT;                              // Array with MCM ADC values (filtered) for tracklet (12bits) //***NEW***
+         Int_t    *fPosLUT;                            // position lookup table **new**
+         UInt_t   *fMCMT;                              // tracklet word for one mcm/trap-chip **new**
+         Int_t   **fZSM;                               // Zero suppression map
           Int_t    *fZSM1Dim;                           // Zero suppression map (1 dimensional projection)
 
+       
+       
+
+
          // Parameter classes
          AliTRDfeeParam *fFeeParam;                    // FEE parameters
          AliTRDSimParam *fSimParam;                    // Simulation parameters
@@ -80,6 +97,7 @@ class AliTRDmcmSim : public TObject {
          void      FilterSimDeConvExpMI(Int_t *source, Double_t *target, Int_t n);
          void      FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n);
          void      FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout, Double_t norm, Double_t lambda, Int_t n);
+         void      FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp);
 
   ClassDef(AliTRDmcmSim,3)
 
index 16e1a7b..08421c3 100644 (file)
@@ -153,6 +153,8 @@ Bool_t AliTRDrawData::Digits2Raw(AliTRDdigitsManager *digitsManager)
 
   // Buffer to temporary store half chamber data
   UInt_t     *hcBuffer    = new UInt_t[kMaxHcWords];
+  
+  Bool_t newEvent = kFALSE;  // only for correct readout tree
 
   // sect is same as iDDL, so I use only sect here.
   for (Int_t sect = 0; sect < fGeo->Nsect(); sect++) { 
@@ -170,11 +172,14 @@ Bool_t AliTRDrawData::Digits2Raw(AliTRDdigitsManager *digitsManager)
     // Reset payload byte size (payload does not include header).
     Int_t npayloadbyte = 0;
 
+    
+
     // GTU common data header (5x4 bytes per super module, shows link mask)
     for( Int_t cham = 0; cham < fGeo->Ncham(); cham++ ) {
       UInt_t gtuCdh = (UInt_t)(0xe << 28);
       for( Int_t plan = 0; plan < fGeo->Nplan(); plan++) {
        Int_t iDet = fGeo->GetDetector(plan, cham, sect);
+       
        // If chamber status is ok, we assume that the optical link is also OK.
         // This is shown in the GTU link mask.
        if ( AliTRDcalibDB::Instance()->GetChamberStatus(iDet) )
@@ -189,10 +194,10 @@ Bool_t AliTRDrawData::Digits2Raw(AliTRDdigitsManager *digitsManager)
       for( Int_t plan = 0; plan < fGeo->Nplan(); plan++) {
 
         Int_t iDet = fGeo->GetDetector(plan,cham,sect);
-
-        // Get the digits array
+       if (iDet == 0) newEvent = kTRUE; // it is expected that each event has at least one tracklet; this is only needed for correct readout tree
+       // Get the digits array
         AliTRDdataArrayS *digits = (AliTRDdataArrayS *) digitsManager->GetDigits(iDet);
-        if (digits->HasData()) {
+        if (digits->HasData() ) {  // second part is new!! and is for indicating a new event
 
           digits->Expand();
 
@@ -204,7 +209,10 @@ Bool_t AliTRDrawData::Digits2Raw(AliTRDdigitsManager *digitsManager)
             hcwords = ProduceHcDataV1andV2(digits,0,iDet,hcBuffer,kMaxHcWords);
          }
          if ( rv == 3 ) { 
-            hcwords = ProduceHcDataV3     (digits,0,iDet,hcBuffer,kMaxHcWords);
+   
+             // hcwords = ProduceHcDataV3     (digits,0,iDet,hcBuffer,kMaxHcWords,newEvent);
+             hcwords = ProduceHcDataV3     (digits,0,iDet,hcBuffer,kMaxHcWords);
+           if(newEvent == kTRUE) newEvent = kFALSE;
          }
 
           of->WriteBuffer((char *) hcBuffer, hcwords*4);
@@ -215,7 +223,9 @@ Bool_t AliTRDrawData::Digits2Raw(AliTRDdigitsManager *digitsManager)
             hcwords = ProduceHcDataV1andV2(digits,1,iDet,hcBuffer,kMaxHcWords);
          }
          if ( rv >= 3 ) {
-            hcwords = ProduceHcDataV3     (digits,1,iDet,hcBuffer,kMaxHcWords);
+          
+             //hcwords = ProduceHcDataV3     (digits,1,iDet,hcBuffer,kMaxHcWords,newEvent);
+             hcwords = ProduceHcDataV3     (digits,1,iDet,hcBuffer,kMaxHcWords);
          }
 
           of->WriteBuffer((char *) hcBuffer, hcwords*4);
@@ -426,8 +436,8 @@ Int_t AliTRDrawData::ProduceHcDataV1andV2(AliTRDdataArrayS *digits, Int_t side
 }
 
 //_____________________________________________________________________________
-Int_t AliTRDrawData::ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side
-                                  , Int_t det, UInt_t *buf, Int_t maxSize)
+//Int_t AliTRDrawData::ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side , Int_t det, UInt_t *buf, Int_t maxSize, Bool_t newEvent)
+Int_t AliTRDrawData::ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side , Int_t det, UInt_t *buf, Int_t maxSize)
 {
   //
   // This function simulates: Raw Version == 3 (Zero Suppression Prototype)
@@ -443,8 +453,10 @@ Int_t AliTRDrawData::ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side
   const Int_t kNTBin = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
   Int_t       kCtype = 0;                       // Chamber type (0:C0, 1:C1)
   //Int_t          iEv = 0xA;                     // Event ID. Now fixed to 10, how do I get event id?
-  UInt_t           x = 0;                       // General used number
-  Int_t           rv = fFee->GetRAWversion();
+
+
+  Bool_t tracklet_on = fFee->GetTracklet();     // **new**
 
   // Check the nCol and nRow.
   if ((nCol == 144) && 
@@ -460,95 +472,115 @@ Int_t AliTRDrawData::ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side
   AliDebug(1,Form("Producing raw data for sect=%d plan=%d cham=%d side=%d"
                  ,sect,plan,cham,side));
 
-  // Tracklet should be processed here but not implemented yet
-
-  // Write end of tracklet marker
-  if (nw < maxSize) {
-    buf[nw++] = kEndoftrackletmarker;
-  } 
-  else {
-    of++;
-  }
+  AliTRDmcmSim** mcm = new AliTRDmcmSim*[(kCtype + 3)*(fGeo->MCMmax())];
 
-  // Half Chamber header
-  // h[0] (there are 3 HC header)
-  Int_t minorv = 0;    // The minor version number
-  Int_t add    = 2;    // The number of additional header words to follow
-  x = (1<<31) | (rv<<24) | (minorv<<17) | (add<<14) | (sect<<9) | (plan<<6) | (cham<<3) | (side<<2) | 1;
-  if (nw < maxSize) {
-    buf[nw++] = x; 
-  }
-  else {
-    of++;
-  }
-  // h[1]
-  Int_t bcCtr   = 99; // bunch crossing counter. Here it is set to 99 always for no reason
-  Int_t ptCtr   = 15; // pretrigger counter. Here it is set to 15 always for no reason
-  Int_t ptPhase = 11; // pretrigger phase. Here it is set to 11 always for no reason
-  x = (bcCtr<<16) | (ptCtr<<12) | (ptPhase<<8) | ((kNTBin-1)<<2) | 1;
-  if (nw < maxSize) {
-    buf[nw++] = x; 
-  }
-  else {
-    of++;
-  }
-  // h[2]
-  Int_t pedSetup       = 1;  // Pedestal filter setup (0:1). Here it is always 1 for no reason
-  Int_t gainSetup      = 1;  // Gain filter setup (0:1). Here it is always 1 for no reason
-  Int_t tailSetup      = 1;  // Tail filter setup (0:1). Here it is always 1 for no reason
-  Int_t xtSetup        = 0;  // Cross talk filter setup (0:1). Here it is always 0 for no reason
-  Int_t nonlinSetup    = 0;  // Nonlinearity filter setup (0:1). Here it is always 0 for no reason
-  Int_t bypassSetup    = 0;  // Filter bypass (for raw data) setup (0:1). Here it is always 0 for no reason
-  Int_t commonAdditive = 10; // Digital filter common additive (0:63). Here it is always 10 for no reason
-  x = (pedSetup<<31) | (gainSetup<<30) | (tailSetup<<29) | (xtSetup<<28) | (nonlinSetup<<27)
-    | (bypassSetup<<26) | (commonAdditive<<20) | 1;
-  if (nw < maxSize) {
-    buf[nw++] = x; 
-  }
-  else {
-    of++;
+  // in case no tracklet-words are processed: write the tracklet-endmarker as well as all additional words immediately and write 
+  // raw-data in one go; if tracklet-processing is enabled, first all tracklet-words of a half-chamber have to be processed before the
+  // additional words (tracklet-endmarker,headers,...)are written. Raw-data is written in a second loop;
+  
+  if (!tracklet_on) {
+      WriteIntermediateWords(buf,nw,of,maxSize,det,side); 
   }
-
+  
   // Scan for ROB and MCM
-  for (Int_t iRobRow = 0; iRobRow < (kCtype + 3); iRobRow++ ) {
+  // scanning direction such, that tracklet-words are sorted in ascending z and then in ascending y order
+  // ROB numbering on chamber and MCM numbering on ROB increase with decreasing z and increasing y
+  for (Int_t iRobRow =  (kCtype + 3)-1; iRobRow >= 0; iRobRow-- ) {
     Int_t iRob = iRobRow * 2 + side;
-    for (Int_t iMcm = 0; iMcm < fGeo->MCMmax(); iMcm++ ) {
-
-      AliTRDmcmSim *mcm = new AliTRDmcmSim();
-      mcm->Init( det, iRob, iMcm );
-      Int_t padrow = mcm->GetRow();
+    // MCM on one ROB
+    for (Int_t iMcmRB = 0; iMcmRB < fGeo->MCMmax(); iMcmRB++ ) {
+       Int_t iMcm = 16 - 4*(iMcmRB/4 + 1) + (iMcmRB%4);
+       Int_t entry = iRobRow*(fGeo->MCMmax()) + iMcm;
+       
+       mcm[entry] = new AliTRDmcmSim();
+       //mcm[entry]->Init( det, iRob, iMcm , newEvent);
+       mcm[entry]->Init( det, iRob, iMcm);
+       //if (newEvent == kTRUE) newEvent = kFALSE; // only one mcm is concerned with new event
+       Int_t padrow = mcm[entry]->GetRow();
+
+       // Copy ADC data to MCM simulator
+       for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
+           Int_t padcol = mcm[entry]->GetCol( iAdc );
+           if ((padcol >=    0) && (padcol <  nCol)) {
+               for (Int_t iT = 0; iT < kNTBin; iT++) { 
+                   mcm[entry]->SetData( iAdc, iT, digits->GetDataUnchecked( padrow, padcol, iT) );
+               } 
+           } 
+           else {  // this means it is out of chamber, and masked ADC
+               mcm[entry]->SetDataPedestal( iAdc );
+           }
+       }
 
-      // Copy ADC data to MCM simulator
-      for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
-       Int_t padcol = mcm->GetCol( iAdc );
-       if ((padcol >=    0) && (padcol <  nCol)) {
-         for (Int_t iT = 0; iT < kNTBin; iT++) { 
-           mcm->SetData( iAdc, iT, digits->GetDataUnchecked( padrow, padcol, iT) );
-         } 
-       } else {  // this means it is out of chamber, and masked ADC
-         mcm->SetDataPedestal( iAdc );
+       // Simulate process in MCM
+       mcm[entry]->Filter();     // Apply filter
+       mcm[entry]->ZSMapping();  // Calculate zero suppression mapping
+
+       if (tracklet_on) {
+           mcm[entry]->Tracklet(); 
+           Int_t tempNw =  mcm[entry]->ProduceTrackletStream( &buf[nw], maxSize - nw );
+           //Int_t tempNw = 0;
+           if( tempNw < 0 ) {
+               of += tempNw;
+               nw += maxSize - nw;
+               AliError(Form("Buffer overflow detected. Please increase the buffer size and recompile."));
+           } else {
+               nw += tempNw;
+           }
        }
-      }
-      // Simulate process in MCM
-      mcm->Filter();     // Apply filter
-      mcm->ZSMapping();  // Calculate zero suppression mapping
-      //mcm->DumpData( "trdmcmdata.txt", "RFZS" ); // debugging purpose
-
-      // Write MCM data to buffer
-      Int_t tempNw =  mcm->ProduceRawStream( &buf[nw], maxSize - nw );
-      if( tempNw < 0 ) {
-       of += tempNw;
-       nw += maxSize - nw;
-       AliError(Form("Buffer overflow detected. Please increase the buffer size and recompile."));
-      } else {
-       nw += tempNw;
-      }
+       // no tracklets: write raw-data already in this loop
+       else {
+            // Write MCM data to buffer
+           Int_t tempNw =  mcm[entry]->ProduceRawStream( &buf[nw], maxSize - nw );
+           if( tempNw < 0 ) {
+               of += tempNw;
+               nw += maxSize - nw;
+               AliError(Form("Buffer overflow detected. Please increase the buffer size and recompile."));
+           } else {
+               nw += tempNw;
+           }
+           
+           delete mcm[entry];
+       }
+           
+       
 
-      delete mcm;
 
+       //mcm->DumpData( "trdmcmdata.txt", "RFZS" ); // debugging purpose
     }
   }
 
+  // if tracklets are switched on, raw-data can be written only after all tracklets
+  if (tracklet_on) {
+      WriteIntermediateWords(buf,nw,of,maxSize,det,side); 
+  
+  
+      // Scan for ROB and MCM
+      for (Int_t iRobRow =  (kCtype + 3)-1; iRobRow >= 0; iRobRow-- ) {
+         //Int_t iRob = iRobRow * 2 + side;
+         // MCM on one ROB
+         for (Int_t iMcmRB = 0; iMcmRB < fGeo->MCMmax(); iMcmRB++ ) {
+             Int_t iMcm = 16 - 4*(iMcmRB/4 + 1) + (iMcmRB%4);
+             
+             Int_t entry = iRobRow*(fGeo->MCMmax()) + iMcm; 
+             
+             // Write MCM data to buffer
+             Int_t tempNw =  mcm[entry]->ProduceRawStream( &buf[nw], maxSize - nw );
+             if( tempNw < 0 ) {
+                 of += tempNw;
+                 nw += maxSize - nw;
+                 AliError(Form("Buffer overflow detected. Please increase the buffer size and recompile."));
+             } else {
+                 nw += tempNw;
+             }
+             
+             delete mcm[entry];
+         
+         }
+      }
+  }
+
+  delete [] mcm;
+  
   // Write end of raw data marker
   if (nw < maxSize) {
     buf[nw++] = kEndofrawdatamarker; 
@@ -560,6 +592,7 @@ Int_t AliTRDrawData::ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side
     AliError("Buffer overflow. Data is truncated. Please increase buffer size and recompile.");
   }
 
+
   return nw;
 
 }
@@ -618,6 +651,66 @@ AliTRDdigitsManager *AliTRDrawData::Raw2Digits(AliRawReader *rawReader)
   return digitsManager;
 }
 
+
+//_____________________________________________________________________________
+void AliTRDrawData::WriteIntermediateWords(UInt_t* buf, Int_t& nw, Int_t& of, const Int_t& maxSize, const Int_t& det, const Int_t& side) {
+    
+    Int_t         plan = fGeo->GetPlane( det );   // Plane
+    Int_t         cham = fGeo->GetChamber( det ); // Chamber
+    Int_t         sect = fGeo->GetSector( det );  // Sector (=iDDL)
+    Int_t           rv = fFee->GetRAWversion();
+    const Int_t kNTBin = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+    UInt_t           x = 0;
+
+    // Write end of tracklet marker
+    if (nw < maxSize) {
+       buf[nw++] = kEndoftrackletmarker;
+    } 
+    else {
+       of++;
+    }
+    
+  // Half Chamber header
+  // h[0] (there are 3 HC header)
+    Int_t minorv = 0;    // The minor version number
+    Int_t add    = 2;    // The number of additional header words to follow
+    x = (1<<31) | (rv<<24) | (minorv<<17) | (add<<14) | (sect<<9) | (plan<<6) | (cham<<3) | (side<<2) | 1;
+    if (nw < maxSize) {
+       buf[nw++] = x; 
+    }
+    else {
+       of++;
+    }
+    // h[1]
+    Int_t bcCtr   = 99; // bunch crossing counter. Here it is set to 99 always for no reason
+    Int_t ptCtr   = 15; // pretrigger counter. Here it is set to 15 always for no reason
+    Int_t ptPhase = 11; // pretrigger phase. Here it is set to 11 always for no reason
+    x = (bcCtr<<16) | (ptCtr<<12) | (ptPhase<<8) | ((kNTBin-1)<<2) | 1;
+    if (nw < maxSize) {
+       buf[nw++] = x; 
+    }
+    else {
+       of++;
+    }
+    // h[2]
+    Int_t pedSetup       = 1;  // Pedestal filter setup (0:1). Here it is always 1 for no reason
+    Int_t gainSetup      = 1;  // Gain filter setup (0:1). Here it is always 1 for no reason
+    Int_t tailSetup      = 1;  // Tail filter setup (0:1). Here it is always 1 for no reason
+    Int_t xtSetup        = 0;  // Cross talk filter setup (0:1). Here it is always 0 for no reason
+    Int_t nonlinSetup    = 0;  // Nonlinearity filter setup (0:1). Here it is always 0 for no reason
+    Int_t bypassSetup    = 0;  // Filter bypass (for raw data) setup (0:1). Here it is always 0 for no reason
+    Int_t commonAdditive = 10; // Digital filter common additive (0:63). Here it is always 10 for no reason
+    x = (pedSetup<<31) | (gainSetup<<30) | (tailSetup<<29) | (xtSetup<<28) | (nonlinSetup<<27)
+       | (bypassSetup<<26) | (commonAdditive<<20) | 1;
+    if (nw < maxSize) {
+       buf[nw++] = x; 
+    }
+    else {
+       of++;
+    } 
+}
+
+
 //_____________________________________________________________________________
 AliTRDdigitsManager *AliTRDrawData::Raw2DigitsOLD(AliRawReader *rawReader)
 {
@@ -714,3 +807,7 @@ AliTRDdigitsManager *AliTRDrawData::Raw2DigitsOLD(AliRawReader *rawReader)
   return digitsManager;
 
 }
+
+
+
+
index f91004b..2335f94 100644 (file)
@@ -41,7 +41,9 @@ class AliTRDrawData : public TObject {
 
   virtual Bool_t       Digits2Raw(AliTRDdigitsManager* digitsManager); // for fRawVersion > 0
   virtual Int_t        ProduceHcDataV1andV2(AliTRDdataArrayS *digits, Int_t side, Int_t det, UInt_t *buf, Int_t maxSize);
+  // virtual Int_t        ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side, Int_t det, UInt_t *buf, Int_t maxSize, Bool_t newEvent);
   virtual Int_t        ProduceHcDataV3(AliTRDdataArrayS *digits, Int_t side, Int_t det, UInt_t *buf, Int_t maxSize);
+          void         WriteIntermediateWords(UInt_t *buf, Int_t& nw, Int_t& of, const Int_t& maxSize, const Int_t& det, const Int_t& side); // writes tracklet-endmarker and additional words between tracklet and raw-data
   
   AliTRDgeometry      *fGeo;            //! Geometry
   AliTRDfeeParam      *fFee;            //! Fee Parameters
diff --git a/TRD/AliTRDtrapAlu.cxx b/TRD/AliTRDtrapAlu.cxx
new file mode 100644 (file)
index 0000000..b10d0ab
--- /dev/null
@@ -0,0 +1,454 @@
+#include "AliTRDtrapAlu.h"\r
+\r
+ClassImp(AliTRDtrapAlu)\r
+\r
+//usage of the class: \r
+//declaration of class instances:   AliTRDtrapAlu a,b,c;\r
+//initialization:                  a.Init(2,11); b.Init(4,4); c.Init(5,4);\r
+//assigning values:                a.AssignDouble(5.7); b.AssignInt(3);(you can also do b.AssignDouble(3) with same effect);\r
+//calculation:                     c = a*b;  \r
+//test if c has right value:       c.WriteWord();\r
+//don't declare pointers; operators not overridden for pointer types; you have to dereference yourself;\r
+//use operators +,-,*,/ only with instances of the class; don't do things like c=a*2 but rather b.AssignInt(2); c=a*b;\r
+\r
+\r
+\r
+  \r
+AliTRDtrapAlu::AliTRDtrapAlu():TObject()\r
+\r
+  ,fValue(0)\r
+  ,fPreCom(0)\r
+  ,fPostCom(0)\r
+  ,fuRestriction(0)\r
+  ,flRestriction(0)\r
+  ,fSigned(kFALSE)\r
+\r
+{\r
+  \r
+ // default constructor\r
\r
\r
+}\r
+\r
+\r
+AliTRDtrapAlu::~AliTRDtrapAlu(){\r
+  //destructor\r
+}\r
+\r
+\r
+\r
+void AliTRDtrapAlu::Init(const Int_t& precom, const Int_t& postcom, const Int_t& lRestriction, const Int_t& uRestriction){\r
+   // initialization: characterizes the bit-word (nr of pre- and post-comma bits, boundaries)\r
+   fPostCom = postcom;\r
+   fPreCom  = precom;\r
+   fValue   = 0;       //currently, re-initialization kills the value\r
+   fSigned  = kFALSE;\r
+\r
+   if (fPreCom + fPostCom > 32 || fPreCom > 31) {fPreCom = 1; fPostCom = 0;return;} // prevent pre-comma part exceeding 31 spaces\r
+   if (fPreCom  <= 0) {fPreCom  = 1;}\r
+   if (fPostCom <  0) {fPostCom = 0;}\r
+   \r
+   Int_t lut = LUT(fPreCom + fPostCom)-1;\r
+   if (uRestriction <= -1 || uRestriction > lut) {fuRestriction = lut;}\r
+   else {fuRestriction = uRestriction;}\r
+   if (lRestriction <= -1 || lRestriction > fuRestriction) {flRestriction = -lut;}\r
+   else {flRestriction = lRestriction;}\r
+   // up to now you can only choose a non-negative lower restriction (e.g. if you want your values to be >=0) ; can't deal with asymmetric borders; have to be implemented if needed\r
+}\r
+\r
+\r
+Double_t AliTRDtrapAlu::GetValueWhole() const { \r
+   // get the actual value (respecting pre- and post-comma parts) in integer-description\r
+   Double_t valPre = (Double_t)(fValue>>fPostCom);\r
+   Double_t valPost = 0.0;\r
+   for(Int_t i = 0; i<=fPostCom-1; i++){\r
+     Double_t num = (fValue>>i)&1;\r
+     Double_t denom = LUT(fPostCom-i);\r
+     valPost = valPost + num/denom;\r
+   }\r
+   Double_t val = valPre + valPost;\r
+   return val;\r
+ }\r
+\r
+\r
+void AliTRDtrapAlu::WriteWord(){\r
+  // for debugging purposes\r
+  printf("bit-word: ");\r
+  if (fSigned == true) printf("-");\r
+  for(Int_t i = fPostCom + fPreCom - 1; i >= fPostCom; i--){  //read from behind in order to write the word from left to right\r
+    printf("%d",(fValue>>i) & 1);\r
+  }\r
+  printf(".");\r
+  for (Int_t j = fPostCom - 1; j >= 0; j--){\r
+    printf("%d",(fValue>>j) & 1);\r
+  }\r
+  printf("\n");\r
+         \r
+}\r
+\r
+\r
+\r
+\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::AssignInt(const Int_t& first){  \r
+  // assign an integer\r
+\r
+  // parameter "first" is an integer for the pre-comma part (not UInt in order to match the error case first<0)\r
+  fSigned = kFALSE;\r
+  Int_t exponent = fPreCom + fPostCom;\r
+\r
+    \r
+  if (first<0) {\r
+    fValue  = 0;                      //setting fValue to 0; first should not be negative\r
+    fValue  = fValue & 0;\r
+    return *this;\r
+  }\r
+\r
+  if (CheckUSize(first<<fPostCom) == kFALSE){\r
+    \r
+    //setting fValue to maximum; first was to big\r
+    fValue  = fuRestriction;\r
+    fValue  = fValue & LUT(exponent)-1;\r
+    return *this;\r
+  }\r
+\r
+  if (CheckLSize(first<<fPostCom) == kFALSE){\r
+    \r
+    //setting fValue to minimum; first was to small\r
+    fValue  = flRestriction;\r
+    fValue  = fValue & LUT(exponent)-1;\r
+    return *this;\r
+  }\r
+\r
+  \r
+  fValue  = first;\r
+  fValue  = fValue<<fPostCom; \r
+  fValue  = fValue & LUT(exponent)-1;\r
\r
+  return *this;\r
+    \r
+}\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::AssignDouble(const  Double_t& first){\r
+  // assign a double\r
\r
+  fSigned = kFALSE;\r
+  Int_t exponent           = fPreCom + fPostCom;\r
+  Int_t firstPre          = 0;  //integer part of first\r
+  Int_t firstPost         = 0;  //comma part of first (cut off with enough accuracy\r
+  Int_t c                  = 0;\r
+  Double_t firstPreFloat = 0;\r
+  \r
+  \r
+  Int_t power1 = LUT(exponent);\r
+  \r
+  firstPre       = (Int_t)first;\r
+  firstPreFloat = firstPre;\r
+  \r
+  if(firstPre < 0){\r
+    fValue  = 0;\r
+    fValue = fValue & 0;\r
+    return *this;\r
+  }\r
+  \r
+  if(CheckUSize((Int_t)(first*LUT(fPostCom))) == kFALSE){\r
+    \r
+    //fValue  = MakePower(2,fPreCom) - 1;\r
+    fValue  = fuRestriction;\r
+    fValue  = fValue & (power1 - 1);\r
+    return *this;\r
+  }\r
+  \r
+  if(CheckLSize((Int_t)(first*LUT(fPostCom))) == kFALSE){\r
+    \r
+    //fValue  = MakePower(2,fPreCom) - 1;\r
+    fValue  = flRestriction;\r
+    fValue  = fValue & (power1 - 1);\r
+    return *this;\r
+  }\r
+  \r
+\r
+  fValue = firstPre;\r
+  \r
+  //get post comma part with adequate accuracy\r
+  firstPost = (Int_t)((first - firstPreFloat)*LUT(fPostCom));\r
+  for(Int_t i = 1; i <= fPostCom; i++) {\r
+    c = (firstPost>>(fPostCom - i)) & 1;\r
+    fValue  = fValue<<1;\r
+    fValue  = fValue | c;\r
+  }\r
+\r
+  fValue = fValue & (power1 - 1);\r
+  return *this;\r
+}\r
+\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::operator=(const AliTRDtrapAlu& binary){\r
+  // assign an object of type AliTRDtrapAlu\r
+\r
+  Int_t c    = 0;\r
+  //Int_t exponent = fPreCom + fPostCom;\r
+  \r
+  \r
+  Int_t power1 = LUT(fPreCom + fPostCom);\r
+\r
+  fValue          = binary.GetValue();         // in case this==&binary : binary's values are overwritten\r
+  Int_t diffPost = binary.GetPost()-fPostCom;\r
+  Int_t check     = 0;\r
+  if(diffPost<0) check = fValue<<(-diffPost);\r
+  else check = fValue>>(diffPost);\r
+  if (CheckUSize(check)==kFALSE){    //checking size of pre-comma part\r
+    \r
+    //setting fValue to maximum\r
+      \r
+           \r
+    fValue  = fuRestriction;         // fuRestriction >= 0 \r
+    fValue  = fValue & (power1 - 1);\r
+    fSigned = kFALSE;\r
+    return *this;\r
+  }\r
+\r
+  Int_t val = (binary.GetSign()==kFALSE) ? check : -check; \r
+  if (CheckLSize(val)==kFALSE){    //checking size of pre-comma part\r
+    \r
+    //setting fValue to minimum\r
+      \r
+           \r
+    if (flRestriction < 0) {\r
+      fValue  = -flRestriction;\r
+      fSigned = kTRUE;\r
+    }\r
+    else {\r
+      fValue  = flRestriction;\r
+      fSigned = kFALSE;\r
+    }\r
+    fValue  = fValue & (power1 - 1);\r
+    return *this;\r
+  }\r
+  \r
+  if (this == & binary) return *this;\r
+  \r
+  fSigned = kFALSE;\r
+  Int_t iValue = fValue;\r
+  fValue = fValue>>(binary.GetPost());           //only keep the valid pre-comma bits\r
+  \r
+  //append existing post-comma bits to fValue; cut off or add 0 if post-comma numbers don`t match\r
+  for(Int_t i = 1; i <= fPostCom; i++){\r
+    if(i <= (binary.GetPost())){\r
+      c = ((iValue)>>(binary.GetPost()-i)) & 1;\r
+    }\r
+    else{\r
+      c = 0;\r
+    }\r
+    fValue  = fValue<<1;\r
+    fValue  = fValue | c;\r
+  }\r
+  \r
+  fValue = fValue & (power1 - 1);\r
+  fSigned = binary.GetSign();\r
+  return *this;\r
+}\r
+\r
+\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::operator+(AliTRDtrapAlu& binary){ \r
+  // + operator\r
+\r
+  //no const parameter because referenced object will be changed\r
+     \r
+  \r
+  Int_t binPre     = binary.GetPre();\r
+  Int_t binPost    = binary.GetPost();\r
+  Int_t binVal     = binary.GetValue();\r
+  \r
+  Int_t min         = Min(binPost,fPostCom);\r
+  Int_t max         = Max(binPre,fPreCom);\r
+  \r
+  Int_t shift       = binPost - min;\r
+  Int_t add1        = (binVal)>>(shift);    //for addition: cut off at minimum accuracy\r
+  shift             = fPostCom - min;\r
+  Int_t add2        = fValue>>(shift);\r
+  if(binary.GetSign() == kTRUE) add1 = -add1;\r
+  if(fSigned == kTRUE) add2 = -add2;\r
+  Int_t add = add1 + add2;\r
+  \r
+  \r
+  //because the parameter "binary" could be a reference to the object to which Mem() is a reference, do not change Mem() until you have extracted all information from "binary"; otherwise you change the information you would like to read\r
+  Mem().Init(max + 1,min);      //buffer: enough space for pre-comma,post-comma according to accuracy\r
+  Mem().AssignFormatted(Max(add,-add));\r
+  Mem().SetSign(add);\r
+  \r
+\r
+  //Mem().FastInit(max+1,min,add);\r
+  return Mem();\r
+}\r
+\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::operator-(AliTRDtrapAlu& binary){\r
+  // - operator    \r
+\r
+  Int_t binPre    = binary.GetPre();\r
+  Int_t binPost   = binary.GetPost();\r
+  Int_t binVal    = binary.GetValue();\r
+\r
+\r
+  Int_t min      = Min(binPost,fPostCom);\r
+  Int_t max      = Max(binPre,fPreCom);\r
+\r
+  Int_t shift    = binPost - min;\r
+  Int_t sub1 = (binVal)>>(shift); //for addition: cut off at minimum accuracy\r
+  shift = fPostCom - min;\r
+  Int_t sub2 = fValue>>(shift);\r
+  if(binary.GetSign() == kTRUE) sub1 = -sub1;\r
+  if(fSigned  == kTRUE) sub2 = -sub2;\r
+  Int_t sub = sub2 - sub1;     // order of subtraction is important\r
\r
+\r
+  Mem().Init(max + 1,min);      //buffer: enough space for pre-comma, post-comma according to accuracy\r
+  Mem().AssignFormatted(Max(sub,-sub)); \r
+  Mem().SetSign(sub);\r
+  //Mem().FastInit(max+1,min,sub);\r
+  return Mem();\r
+} \r
+\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::operator*(AliTRDtrapAlu& binary){\r
+  // * operator\r
+  \r
+\r
+  Int_t binPre   = binary.GetPre();\r
+  Int_t binPost  = binary.GetPost();\r
+\r
+\r
+  Int_t min      = Min(binPost,fPostCom);\r
+  Int_t max      = Max(binPre,fPreCom);\r
+\r
+  \r
+  Int_t mult1 = binary.GetValue();\r
+  Int_t mult2 = fValue;\r
+  Int_t shift  = (Int_t)(fPostCom + binPost - min);\r
+  Double_t fmult1 = (Double_t)mult1;\r
+  Double_t fmult2 = (Double_t)mult2;\r
+  (fmult1 > fmult2) ? fmult1 = fmult1/LUT(shift) : fmult2 = fmult2/LUT(shift);\r
+  \r
+    \r
+  if (binary.GetSign() == kTRUE) fmult1 = -fmult1;\r
+  if (fSigned  == kTRUE) fmult2 = -fmult2;\r
+  Double_t fmult  = fmult1*fmult2;\r
+  Int_t mult = (Int_t)fmult;\r
+  Int_t sign = 1;\r
+  if(mult<0) sign = -1;\r
+  mult = Max(mult,-mult);\r
+  //Int_t shift = fPostCom + binPost - min;\r
+  //mult = mult>>(shift);\r
+  Mem().Init(2 * max + 1, min); // +1 to consider the borrow from the past-comma part; accuracy of past-comma part is determined by the minimum; therefore, for the result not more accuracy is guaranteed\r
+  // be aware that this only works if 2*max+1+min <= 32!! adjusting the pre-comma place to the value would consume too much time\r
+\r
+  Mem().AssignFormatted(mult);\r
+  Mem().SetSign(sign);\r
+  //mult = sign*mult;\r
+  //Mem().FastInit(2*max+1,min,mult);\r
+  return Mem();\r
+}\r
+  \r
+\r
+\r
+AliTRDtrapAlu& AliTRDtrapAlu::operator/(AliTRDtrapAlu& binary){\r
+  // / operator\r
\r
+\r
+  Int_t binPre  = binary.GetPre();\r
+  Int_t binPost = binary.GetPost();\r
+  Int_t min              = Min(binPost,fPostCom);\r
+  Int_t max              = Max(binPre,fPreCom);\r
+  \r
+  Int_t div1             = binary.GetValue(); //value in integer format\r
+  Int_t div2             = fValue;\r
+  \r
+  // this approach does not always work because it can exceed the range of integers\r
+  //Int_t numerator     = div2 * LUT(min);\r
+  Int_t numerator     = div2;\r
+  if (fSigned == kTRUE) numerator = numerator*(-1);\r
+  Int_t denominator   = div1;\r
+  if (binary.GetSign() == kTRUE) denominator = denominator*(-1);\r
+  Double_t fdiv       = 0.0;\r
+  Double_t fLUT       = 0.0;\r
+  Int_t div           = 0;\r
+  \r
+\r
+  if (div1 == 0){\r
+    Mem().Init(max + 1,min);\r
+    Mem().AssignFormatted(LUT(max+min+1)-1); // division by 0: set to max value\r
+    //Mem().FastInit(max+1,min,div1);\r
+    return Mem();\r
+  }\r
+      \r
+  fdiv = (Double_t)numerator/denominator;\r
+  \r
+  Int_t shift = fPostCom - binPost;\r
+  \r
+  if(shift>0){\r
+      //denominator = denominator * LUT(shift);\r
+      fLUT = (Double_t)LUT(min)/LUT(shift);\r
+  }\r
+  else {\r
+      if(shift<0) {\r
+         shift = -shift;\r
+         //numerator =  numerator * LUT(shift);\r
+         fLUT = (Double_t)LUT(min)*LUT(shift);\r
+      }\r
+      else {\r
+         fLUT = (Double_t)LUT(min);\r
+      }\r
+  }\r
+\r
+  fdiv = fdiv*fLUT;\r
+  div = (Int_t)fdiv;\r
+  \r
+  Int_t sign = (div>=0) ? 1 : -1;\r
+  div = Max(div,-div);\r
+  \r
+  // chose min as past-comma part because from a division of integers you can't get only an integer\r
+  Mem().Init(max + 1,min); // max+1+min must <= 32!!\r
+  Mem().SetSign(sign);\r
+  Mem().AssignFormatted(div);\r
+  \r
+  return Mem();\r
+  \r
+}\r
+\r
+\r
+\r
+Int_t AliTRDtrapAlu::MakePower(const Int_t& base,const  Int_t& exponent)const{\r
+// calculate "base" to the power of "exponent"\r
+  Int_t result = 1;\r
+  \r
+    for(Int_t i = 1; i <= exponent; i++){\r
+    result = result * base;\r
+  }\r
+  return result;\r
+}\r
+\r
+\r
+\r
+\r
+Int_t AliTRDtrapAlu::LUT(const Int_t& index){   \r
+  // simple look-up table for base=2\r
+  \r
+  static Bool_t fLUT = kFALSE;\r
+  static Int_t gLUT[30];\r
+  if (fLUT == kFALSE) {\r
+    gLUT[0] = 1;\r
+    for(Int_t i = 1; i<30; i++) {\r
+      gLUT[i] = gLUT[i-1] * 2;\r
+    }\r
+  fLUT = kTRUE;\r
+  } \r
+  if (index >=0 && index < 30){\r
+    return gLUT[index];\r
+  }\r
+  else {\r
+    \r
+    return Mem().MakePower(2,index);\r
+  }\r
+}\r
diff --git a/TRD/AliTRDtrapAlu.h b/TRD/AliTRDtrapAlu.h
new file mode 100644 (file)
index 0000000..b4cb7c5
--- /dev/null
@@ -0,0 +1,160 @@
+#ifndef ALITRDTRAPALU_H\r
+#define ALITRDTRAPALU_H\r
+#include <TObject.h>\r
+\r
+class AliTRDtrapAlu:public TObject {\r
+\r
+//class AliTRDtrapAlu {\r
\r
+ public:\r
+   \r
+               AliTRDtrapAlu();\r
+              //AliTRDtrapAlu(AliTRDtrapAlu& bin); //copy constructor\r
+ virtual      ~AliTRDtrapAlu();\r
+\r
+\r
+ void Init(const Int_t& precom=10, const Int_t& postcom=2, const Int_t& lRestriction = -1, const Int_t& uRestriction = -1);\r
+      \r
+\r
+\r
+Int_t   GetValue ()const{ \r
+    // return the value \r
+    return fValue;\r
+ }\r
+\r
+Int_t   GetSignedValue ()const{  \r
+    // return the value with its sign\r
+    if(fSigned == kFALSE) return fValue;\r
+    else return fValue*(-1);\r
+ }\r
+\r
+Int_t   GetValuePre ()const{\r
+    // return value of pre-comma part as integer\r
+    Int_t valPre = fValue>>fPostCom;\r
+    return valPre;   \r
+}\r
+\r
+Double_t GetValueWhole() const;  \r
+   \r
+\r
+Int_t   GetPre()const{\r
+    // return nr of pre-comma bits\r
+    return fPreCom;\r
+}\r
+\r
+Int_t   GetPost()const{\r
+    // return nr of past-comma bits\r
+    return fPostCom;\r
+}\r
+\r
+Bool_t  GetSign()const{\r
+    // return true if signed\r
+    if(fSigned == kTRUE) return kTRUE;\r
+    return kFALSE;\r
+}\r
+\r
+      \r
+Bool_t  CheckUSize(const Int_t& val)const{\r
+    // compare value to the upper restriction\r
+    if(val>fuRestriction) return kFALSE;\r
+    return kTRUE;\r
+}\r
+\r
+Bool_t  CheckLSize(const Int_t& val)const{\r
+    // compare value to the lower restriction\r
+    if(val<flRestriction) return kFALSE;\r
+    return kTRUE;\r
+}\r
+\r
+\r
+void AssignFormatted(const Int_t& formVal){ \r
+// assign a value with proper format; assigns formVal directly to fValue; better not use explicitely\r
+    fValue  = formVal;\r
+    //fValue  = fValue & (LUT(fPreCom + fPostCom) - 1); // no cut-off wanted\r
+} \r
+\r
+void SetSign(const Int_t& s){\r
+    // sets the sign\r
+    if(s >= 0) fSigned = kFALSE;\r
+    if(s <  0) fSigned = kTRUE;\r
+}\r
+\r
+\r
+void   WriteWord();  \r
+\r
+\r
+AliTRDtrapAlu& AssignInt(const  Int_t& first);     // in case a decimal integer is assigned to a binary; \r
+AliTRDtrapAlu& AssignDouble(const  Double_t& first);  // change "Double_t" into "Float_t"\r
+AliTRDtrapAlu& operator=(const AliTRDtrapAlu& binary);\r
+       \r
+AliTRDtrapAlu& operator+(AliTRDtrapAlu& binary); //binary is not const, because in a+(b*c) binary is reference to the object, to which Mem() is also a reference and this object is changed\r
+\r
+AliTRDtrapAlu& operator-(AliTRDtrapAlu& binary);\r
+\r
+AliTRDtrapAlu& operator*(AliTRDtrapAlu& binary);\r
+\r
+AliTRDtrapAlu& operator/(AliTRDtrapAlu& binary);\r
+      \r
+\r
+\r
+\r
+\r
+protected:\r
+\r
+// void FastInit(const Int_t& precom = 10, const Int_t&  postcom = 2, const Int_t& formVal = 0); //meant to combine definition of format with integer-value assignment; not to apply by user  \r
+\r
+\r
+//the following two functions encapsulate global static members; can only be changed by member functions (visibility only inside class)\r
+\r
+\r
\r
\r
+ Int_t  MakePower(const Int_t& base=1,const  Int_t& exponent=1)const;\r
+\r
+\r
+ static AliTRDtrapAlu& Mem() { \r
+// a global instance of the class, which is only defined once\r
+   static AliTRDtrapAlu fAuxiliary;\r
+   return fAuxiliary;\r
+ }\r
+\r
+\r
+ static Int_t LUT(const Int_t& index);\r
+       \r
+\r
+\r
+       \r
+const Int_t&  Min(const Int_t& comp1, const Int_t& comp2)const{\r
+    // return the minimum\r
+    if (comp1 <= comp2) return comp1;\r
+    return comp2;\r
+}\r
+\r
+const Int_t&  Max(const Int_t& comp1, const Int_t& comp2)const{\r
+    // return the maximum\r
+    if (comp1 >= comp2) return comp1;\r
+    return comp2;\r
+}\r
+\r
+\r
+       Int_t    fValue;          // the value in integers\r
+       Int_t    fPreCom;         // number of pre-comma bits\r
+       Int_t    fPostCom;        // number of past-comma bits\r
+       Int_t    fuRestriction;   // the upper restriction for the value\r
+       Int_t    flRestriction;   // the lower restriction for the value\r
+       Bool_t   fSigned;         // signed value? \r
+      \r
+      \r
+       \r
+\r
+\r
+\r
+\r
+ClassDef(AliTRDtrapAlu,1)\r
+};\r
+\r
+\r
+\r
+#endif\r
+\r
+\r