]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDmcmSim.cxx
Update of raw data simulation by WooJin
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
index 292235cd3416a65bd4ea69ed987858accc70bb3e..264efaa6c031e8da57b8465f1cc0265fd6c4d3d5 100644 (file)
@@ -101,19 +101,23 @@ The default raw version is 2.
 #include "AliTRDSimParam.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDcalibDB.h"
-
+#include "AliTRDdigitsManager.h"
+#include "AliTRDarrayADC.h"
 // additional for new tail filter and/or tracklet
 #include "AliTRDtrapAlu.h"
 #include "AliTRDpadPlane.h"
+#include "AliTRDtrackletMCM.h"
 
+#include "AliRun.h"
+#include "AliLoader.h"
 
 ClassImp(AliTRDmcmSim)
 
 //_____________________________________________________________________________
 AliTRDmcmSim::AliTRDmcmSim() :TObject()
   ,fInitialized(kFALSE)
-  ,fNextEvent(-1)    //new
-  ,fMaxTracklets(-1) //new
+  ,fNextEvent(-1)    
+  ,fMaxTracklets(-1) 
   ,fChaId(-1)
   ,fSector(-1)
   ,fStack(-1)
@@ -122,12 +126,12 @@ AliTRDmcmSim::AliTRDmcmSim() :TObject()
   ,fMcmPos(-1)
   ,fNADC(-1)
   ,fNTimeBin(-1)
-  ,fRow(-1)
+  ,fRow (-1)
   ,fADCR(NULL)
   ,fADCF(NULL)
-  ,fADCT(NULL)      //new
-  ,fPosLUT(NULL)    //new
-  ,fMCMT(NULL)      //new
+  ,fADCT(NULL)     
+  ,fPosLUT(NULL)    
+  ,fMCMT(NULL)      
   ,fZSM(NULL)
   ,fZSM1Dim(NULL)
   ,fFeeParam(NULL)
@@ -147,8 +151,8 @@ AliTRDmcmSim::AliTRDmcmSim() :TObject()
 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m) 
   :TObject(m)
   ,fInitialized(kFALSE) 
-  ,fNextEvent(-1)    //new
-  ,fMaxTracklets(-1) //new
+  ,fNextEvent(-1)    
+  ,fMaxTracklets(-1) 
   ,fChaId(-1)
   ,fSector(-1)
   ,fStack(-1)
@@ -160,9 +164,9 @@ AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
   ,fRow(-1)
   ,fADCR(NULL)
   ,fADCF(NULL)
-  ,fADCT(NULL)      //new
-  ,fPosLUT(NULL)    //new
-  ,fMCMT(NULL)      //new
+  ,fADCT(NULL)      
+  ,fPosLUT(NULL)    
+  ,fMCMT(NULL)      
   ,fZSM(NULL)
   ,fZSM1Dim(NULL)
   ,fFeeParam(NULL)
@@ -256,23 +260,24 @@ void AliTRDmcmSim::Copy(TObject &m) const
 }
 
 //_____________________________________________________________________________
-//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)
+
+//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 = kFALSE ) // 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!!**
+   
+  fNextEvent     = 0; 
   fFeeParam      = AliTRDfeeParam::Instance();
   fSimParam      = AliTRDSimParam::Instance();
   fCal           = AliTRDcalibDB::Instance();
   fGeo           = new AliTRDgeometry();
   fChaId         = chaId;
   fSector        = fGeo->GetSector( fChaId );
-  fStack         = fGeo->GetChamber( fChaId );
-  fLayer         = fGeo->GetPlane( fChaId );
+  fStack         = fGeo->GetStack( fChaId );
+  fLayer         = fGeo->GetLayer( fChaId );
   fRobPos        = robPos;
   fMcmPos        = mcmPos;
   fNADC          = fFeeParam->GetNadcMcm();
@@ -283,11 +288,11 @@ void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos ) // only for r
 
  
 
-  //if ((fChaId == 0 && fRobPos == 0 && fMcmPos == 0)) {
-  /* if (newEvent == kTRUE) {
-    fNextEvent = 1;
-    }*/
-   fNextEvent = 0;
+  
+  if (newEvent == kTRUE) {
+      fNextEvent = 1;
+  }
+
 
 
   // Allocate ADC data memory if not yet done
@@ -348,17 +353,17 @@ 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();
+  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 kNlayer = fGeo->Nlayer();
 
-  Float_t  *sPRFsmp   = new Float_t[nBin*kNplan];
+  Float_t  *sPRFsmp   = new Float_t[nBin*kNlayer];
   Double_t *sPRFlayer = new Double_t[nBin];
   
   
-  for(Int_t i = 0; i<nBin*kNplan; i++){
+  for(Int_t i = 0; i<nBin*kNlayer; i++){
     
     //printf("%f\n",fCal->GetSampledPRF()[i]);
     sPRFsmp[i] = fCal->GetSampledPRF()[i]; 
@@ -604,6 +609,94 @@ Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
   if( of != 0 ) return -of; else return nw;
 }
 
+//_____________________________________________________________________________
+Int_t AliTRDmcmSim::ProduceRawStreamV2( UInt_t *buf, Int_t maxSize, UInt_t iEv )
+{
+  //
+  // Produce raw 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;
+  //UInt_t  iEv = 0;
+  Int_t   nw  = 0;  // Number of written words
+  Int_t   of  = 0;  // Number of overflowed words
+  Int_t   rawVer   = fFeeParam->GetRAWversion();
+  Int_t **adc;
+  Int_t   nActiveADC = 0;      // number of activated ADC bits in a word
+
+  if( !CheckInitialized() ) return 0;
+
+  if( fFeeParam->GetRAWstoreRaw() ) {
+    adc = fADCR;
+  } else {
+    adc = fADCF;
+  }
+
+  // Produce MCM header : xrrr mmmm eeee eeee eeee eeee eeee 1100
+  //                      x : 0 before , 1 since 10.2007
+  //                      r : Readout board position (Alice numbering)
+  //                      m : MCM posi
+  //                      e : Event counter from 1
+  //x = (1<<31) | ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
+  x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
+  if (nw < maxSize) {
+    buf[nw++] = x;
+       //printf("\nMCM header: %X ",x);
+  }
+  else {
+    of++;
+  }
+
+  // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
+  //                           n : unused , c : ADC count, m : selected ADCs
+  if( rawVer >= 3 ) {
+    x = 0;
+    for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
+      if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
+               x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
+               nActiveADC++;           // number of 1 in mmm....m
+      }
+    }
+       x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
+       //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
+
+    if (nw < maxSize) {
+      buf[nw++] = x;
+         //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
+    }
+    else {
+      of++;
+    }
+  }
+
+  // Produce ADC data. 3 timebins are packed into one 32 bits word
+  // In this version, different ADC channel will NOT share the same word
+
+  UInt_t aa=0, a1=0, a2=0, a3=0;
+
+  for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
+    if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
+    aa = !(iAdc & 1) + 2;      // 3 for the even ADC channel , 2 for the odd ADC channel
+    for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
+      a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
+      a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
+      a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
+      x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
+      if (nw < maxSize) {
+       buf[nw++] = x;
+       //printf("%08X ",x);
+      }
+      else {
+       of++;
+      }
+    }
+  }
+
+  if( of != 0 ) return -of; else return nw;
+}
+
 //_____________________________________________________________________________
 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
 {
@@ -663,6 +756,8 @@ void AliTRDmcmSim::Filter()
 //_____________________________________________________________________________
 void AliTRDmcmSim::FilterPedestal()
 {
+
+     
   //
   // Apply pedestal filter
   //
@@ -747,8 +842,9 @@ void AliTRDmcmSim::FilterTail()
     break;
   }
 
-  delete dtarg;
-  delete itarg;
+  delete [] dtarg;
+  delete [] itarg;
+
 }
 
 //_____________________________________________________________________________
@@ -1249,20 +1345,20 @@ void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, I
 
 
 
-
 //__________________________________________________________________________________
 
 
 // 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
+//    currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
+// -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with          information about the Tracklet position (MCM, channel number)
 
-// 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
+// The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
 
 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 you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
 
   if(!CheckInitialized()){ return; }
   
@@ -1277,13 +1373,13 @@ void AliTRDmcmSim::Tracklet(){
     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
+       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
+  // the online ordering of mcm's is reverse to the TRAP-manual-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];
@@ -1343,7 +1439,7 @@ void AliTRDmcmSim::Tracklet(){
     delete[] isource;
   }
   
-  //the following values should be in AliTRDfeeParam and have to be read in properly
+  //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
   //naming follows conventions in TRAP-manual
   
   
@@ -1417,10 +1513,13 @@ void AliTRDmcmSim::Tracklet(){
   QT0Alu.Init(15,0);
   AliTRDtrapAlu QT1Alu;
   QT1Alu.Init(16,0);
+
+  AliTRDtrapAlu oneAlu;
+  oneAlu.Init(1,8);
  
   
   AliTRDtrapAlu inverseNAlu;
-  inverseNAlu.Init(1,8);        // replaces the LUT for 1/N
+  inverseNAlu.Init(1,8);        // simulates the LUT for 1/N
   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
   MeanChargeAlu.Init(8,0);
   AliTRDtrapAlu TotalChargeAlu;
@@ -1459,7 +1558,7 @@ void AliTRDmcmSim::Tracklet(){
   
   for(Int_t iT = 0; iT < fNTimeBin; iT++){
     
-    if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no? // !!**enable**!!
+    if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no? 
 
     // reset
     Int_t portChannel[4]   = {-1,-1,-1,-1};   
@@ -1470,7 +1569,7 @@ void AliTRDmcmSim::Tracklet(){
     
     Int_t mark = 0;
     
-    filPed.AssignFormatted(ieffped[iT]);   // no size-checking with AssignFormatted; ieffped>=0
+    filPed.AssignFormatted(ieffped[iT]);   // no size-checking when using AssignFormatted; ieffped>=0
     filPed = filPed;                       // this checks the size
     
     ieffped[iT] = filPed.GetValue();
@@ -1516,7 +1615,7 @@ void AliTRDmcmSim::Tracklet(){
     Int_t j = 1; // selection number
     while(i<fNADC-2 && j<=3){
       i = i + 1;
-      if((mark>>(i-1)) & 1 == 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
        
@@ -1552,7 +1651,7 @@ void AliTRDmcmSim::Tracklet(){
     // read three from left side
     Int_t k = fNADC-2;
     while(k>i && j<=6) {
-      if((mark>>(k-1)) & 1 == 1) {
+      if( ((mark>>(k-1)) & 1) == 1) {
        selection[j].iadc  = fNADC-1-k;
        selection[j].value = qsum[fNADC-1-k]>>6;
        
@@ -1722,7 +1821,7 @@ void AliTRDmcmSim::Tracklet(){
   
   // 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
+  const Int_t cTCT = 8;      // joint number of hits;     8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
   
   Int_t mPair   = 0;         // marker for possible tracklet pairs
   Int_t* hitSum = new Int_t[fNADC-3];
@@ -1759,7 +1858,7 @@ void AliTRDmcmSim::Tracklet(){
   // 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) {
+    if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
       selectPair[selNr].iadc  = adcL;
       selectPair[selNr].value = hitSum[adcL];   
       
@@ -1813,11 +1912,11 @@ void AliTRDmcmSim::Tracklet(){
        rpairChannel[i] = -1;
        break;
       }
-      if(rpairChannel[i] == lpairChannel[j]) {
+      /* if(rpairChannel[i] == lpairChannel[j]) {
        lpairChannel[i] = -1;
        rpairChannel[i] = -1;
        break;
-      }
+       }*/
     }
   }
   
@@ -1838,7 +1937,8 @@ void AliTRDmcmSim::Tracklet(){
   Int_t mMeanCharge[4]; 
   Int_t inverseN = 0;
   Double_t invN = 0;
+  Int_t one = 0;
+
   for (Int_t i = 0; i<4; i++){
     mADC[i] = -1;                        // set to invalid number
     mN[i]   =  0;
@@ -1854,16 +1954,20 @@ void AliTRDmcmSim::Tracklet(){
     mMeanCharge[i] = 0;
   }
   
-
-  YAlu.AssignInt(1);
-  Int_t wpad  = YAlu.GetValue();       // 1 with 8 past-comma bits
+  oneAlu.AssignInt(1);
+  one = oneAlu.GetValue();              // one 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
+                                        // 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;
+    
+    YAlu.AssignInt(N[rpairChannel[i]]);
+    Int_t wpad  = YAlu.GetValue();       // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
+
     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
@@ -1900,6 +2004,7 @@ void AliTRDmcmSim::Tracklet(){
     MeanChargeAlu = TotalChargeAlu;
     mMeanCharge[i] = MeanChargeAlu.GetValue();
     
+    // this check is not necessary; it is just for efficiency reasons
     if (N[madc+1] == 0) {
        mX[i]     =   X[madc];
        mXX[i]    =  XX[madc];
@@ -1911,12 +2016,12 @@ void AliTRDmcmSim::Tracklet(){
        
        mX[i]     =   X[madc] +  X[madc+1];
        XAlu.AssignFormatted(mX[i]);
-       XAlu    = XAlu;
+       XAlu      = XAlu;
        mX[i]     = XAlu.GetValue();
        
        mXX[i]    =  XX[madc] + XX[madc+1];
        XXAlu.AssignFormatted(mXX[i]);
-       XXAlu   = XXAlu;
+       XXAlu     = XXAlu;
        mXX[i]    = XXAlu.GetValue();
  
     
@@ -1932,7 +2037,7 @@ void AliTRDmcmSim::Tracklet(){
        YAlu    = YAlu;
        mY[i]     = YAlu.GetSignedValue();
        
-       mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*wpad;
+       mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*one;    // multiplication by one to maintain the data format
        
        if (mXY[i] < 0) {
            XYAlu.AssignFormatted(-mXY[i]);
@@ -1945,7 +2050,7 @@ void AliTRDmcmSim::Tracklet(){
        XYAlu   = XYAlu;
        mXY[i]    = XYAlu.GetSignedValue();
     
-       mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*wpad + wpad*wpad;
+       mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
        if (mYY[i] < 0) {
            YYAlu.AssignFormatted(-mYY[i]);
            YYAlu.SetSign(-1);
@@ -1965,6 +2070,16 @@ void AliTRDmcmSim::Tracklet(){
   // 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
   
+  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  // !!important note: the offset is calculated from hits in the time bin range between tFS and tFE; it corresponds to the value at the height of the time bin tFS which does NOT need to correspond to the upper side of the drift   !!
+  // !!volume (cathode wire plane). The offset cannot be rescaled as long as it is unknown which is the first time bin that contains hits from the drift region and thus to which distance from the cathode plane tFS corresponds.    !!
+  // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
+  // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
+  // !!The corrections to the offset (e.g. no ExB correction applied as offset is supposed to be on top of drift region; however not at anode wire, so some inclination of drifting clusters due to Lorentz angle exists) are only    !!
+  // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
+  // !!The slope however can be converted to a deflection length between electrode and cathode wire plane as it is clear that the drift region is sampled 20 times                                                                    !!
+  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
   // which formats should be chosen?
   AliTRDtrapAlu denomAlu;
   denomAlu.Init(20,8);       
@@ -2021,7 +2136,7 @@ void AliTRDmcmSim::Tracklet(){
     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
+    // up to here: adc-mapping according to TRAP manual 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; 
@@ -2039,8 +2154,8 @@ void AliTRDmcmSim::Tracklet(){
   
   // 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 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
@@ -2063,12 +2178,17 @@ void AliTRDmcmSim::Tracklet(){
 
       
   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
+  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 dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
+  Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
+  Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
+  Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
+  Int_t    nrOfAmplTimeBins  = 2;                           // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
+  Int_t    nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is  to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
+  if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed 
   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 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);
@@ -2080,9 +2200,9 @@ void AliTRDmcmSim::Tracklet(){
   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;
+  // 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: 
@@ -2118,12 +2238,32 @@ void AliTRDmcmSim::Tracklet(){
     numAlu.AssignInt(nCorrectOffset);
     numAlu.SetSign(1);
   }
-  nCorrectOffset = numAlu.GetSignedValue();         
+  nCorrectOffset = numAlu.GetSignedValue();   
+
+  // the Lorentz correction to the offset
+  Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
+  
+
+  lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
+  
+  if(lorCorrectOffset < 0) {
+      numAlu.AssignDouble(-lorCorrectOffset);
+      numAlu.SetSign(-1);
+  }
+  else{
+      numAlu.AssignDouble(lorCorrectOffset);
+      numAlu.SetSign(1);
+  }
+  
+  Int_t mlorCorrectOffset = numAlu.GetSignedValue();
+  
+  
   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
  
   // calculation of slope-correction
 
   // this is only true for tracks coming (approx.) from primary vertex
+  // everything is evaluated for a tracklet covering the whole drift chamber
   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
@@ -2140,10 +2280,11 @@ void AliTRDmcmSim::Tracklet(){
  
   // 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;
+  // note that the fit was only done in the range tFS to tFE, however this range does not need to cover the whole drift region (neither start nor end of it)
+  // however the tracklets are supposed to be a fit in the drift region thus the linear function is stretched to fit the drift region of 30 mm
+  
   
-  Double_t mCorrectSlope = ((Double_t)(fNTimeBin-tFS))*padWidthI/granularitySlope;  // >= 0.0
+  Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
 
   AliTRDtrapAlu correctAlu;
   correctAlu.Init(20,8);
@@ -2166,8 +2307,8 @@ void AliTRDmcmSim::Tracklet(){
     
     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
+    mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset)); 
+    //mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
     
     if (mOffset[i] < 0) {
       numAlu.AssignFormatted(-mOffset[i]);
@@ -2266,72 +2407,13 @@ void AliTRDmcmSim::Tracklet(){
       //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 
+    //bit-word is 2-complement and therefore without sign
     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);
-      }*/
+    UInt_t shift  = 0;
+    UInt_t shift2 = 0;
+       
+       
 
 
     /*printf("mean charge: %d\n",mMeanCharge[i]);
@@ -2341,64 +2423,82 @@ void AliTRDmcmSim::Tracklet(){
     printf("channel: %d\n",mADC[i]);*/
 
     // electron probability (currently not implemented; the mean charge is just scaled)
+    shift = (UInt_t)mMeanCharge[i];
     for(Int_t iBit = 0; iBit < 8; iBit++) {
       bitWord[j]  = bitWord[j]<<1;
-      bitWord[j] |= (mMeanCharge[i]>>(7-iBit))&1;               
+      bitWord[j] |= (shift>>(7-iBit))&1;               
       //printf("0");
     }
-    
+
     // pad row
+    shift = (UInt_t)fRow;
     for(Int_t iBit = 0; iBit < 4; iBit++) {
       bitWord[j]  = bitWord[j]<<1;
-      bitWord[j] |= (fRow>>(3-iBit))&1;
+      bitWord[j] |= (shift>>(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));
-      }
+       shift = (UInt_t)(-mSlope[i]);
+       // shift2 is 2-complement of shift
+       shift2 = 1;
+       for(Int_t iBit = 1; iBit < 7; iBit++) {
+           shift2  = shift2<<1;
+           shift2 |= (1- (((shift)>>(6-iBit))&1) );
+           //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
+       }
+       shift2 = shift2 + 1;
+       //printf("1");
+       for(Int_t iBit = 0; iBit < 7; iBit++) {
+           bitWord[j]  = bitWord[j]<<1;
+           bitWord[j] |= (shift2>>(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++) {
+       shift = (UInt_t)(mSlope[i]);
        bitWord[j]  = bitWord[j]<<1;
-       bitWord[j] |= (mSlope[i]>>(6-iBit))&1;
-       //printf("%d",(mSlope[i]>>(6-iBit))&1);
-      }
+       bitWord[j]   |= 0;
+       //printf("0");
+       for(Int_t iBit = 1; iBit < 7; iBit++) {
+           bitWord[j]  = bitWord[j]<<1;
+           bitWord[j] |= (shift>>(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");
+       shift = (UInt_t)(-mOffset[i]);
+       shift2 = 1;
        for(Int_t iBit = 1; iBit < 13; iBit++) {
-           bitWord[j]  = bitWord[j]<<1;
-           bitWord[j] |= (1-((-mOffset[i])>>(12-iBit))&1);
+           shift2  = shift2<<1;
+           shift2 |= (1-(((shift)>>(12-iBit))&1));
            //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
        }
+       shift2 = shift2 + 1;
+       //printf("1");
+       for(Int_t iBit = 0; iBit < 13; iBit++) {
+           bitWord[j]  = bitWord[j]<<1;
+           bitWord[j] |= (shift2>>(12-iBit))&1;
+           //printf("%d",(1-((-mSlope[i])>>(6-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);
-      }
+       shift = (UInt_t)mOffset[i];
+       bitWord[j] = bitWord[j]<<1;
+       bitWord[j]   |= 0; 
+       //printf("0");
+       for(Int_t iBit = 1; iBit < 13; iBit++) {
+           bitWord[j]  = bitWord[j]<<1;
+           bitWord[j] |= (shift>>(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];
@@ -2424,255 +2524,28 @@ void AliTRDmcmSim::Tracklet(){
 
   delete padPlane;
 
-
-
-/*
+//if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
+       
+  if (!fFeeParam->GetMCTrackletOutput()) 
+      return;
  
-
-  // 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) {
+  if (wordnr == 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);
-    }
-  }
+  UInt_t       *trackletWord = new UInt_t[fMaxTracklets];
+  Int_t        *adcChannel   = new Int_t[fMaxTracklets];
+  Int_t        *trackRef     = new Int_t[fMaxTracklets];
 
-      
-  
-  // fill the tracklet word; here wordnr > 0
+  Int_t u = 0;
 
-  trackletWord = new UInt_t[fMaxTracklets];
-  adcChannel   = new Int_t[fMaxTracklets];
+  AliTRDdigitsManager *digman = new AliTRDdigitsManager();
+  digman->ReadDigits(AliRunLoader::GetRunLoader()->GetLoader("TRDLoader")->TreeD());
+  digman->SetUseDictionaries(kTRUE);
+  AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
 
   for (Int_t j = 0; j < fMaxTracklets; j++) {
       Int_t i = order[j];
@@ -2681,57 +2554,211 @@ void AliTRDmcmSim::Tracklet(){
       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];
+
+// Finding label of MC track
+         TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
+         Int_t track[3];
+         Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
+         Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
+         Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
+         Int_t det = 30 * fSector + 6 * fStack + fLayer;
+         for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
+             track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
+             track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
+             track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
+             hTrkRef->Fill(track[0]);
+             if (track[1] != track[0] && track[1] != -1)
+                 hTrkRef->Fill(track[1]);
+             if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
+                 hTrkRef->Fill(track[2]);
+             if (padcol_ngb >= 0) {
+                 track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
+                 track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
+                 track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
+                 hTrkRef->Fill(track[0]);
+                 if (track[1] != track[0] && track[1] != -1)
+                     hTrkRef->Fill(track[1]);
+                 if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
+                     hTrkRef->Fill(track[2]);
+             }
+         }
+         trackRef[u] = hTrkRef->GetMaximumBin() - 1;
+         delete hTrkRef;
          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();
+  AliDataLoader *dl = AliRunLoader::GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
+  if (!dl) {
+    AliError("Could not get the tracklets data loader!");
+  }
+  else {
+    TTree *trackletTree = dl->Tree();
+    if (!trackletTree)
+      dl->MakeTree();
+    trackletTree = dl->Tree();
+
+   AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM(); 
+   TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
+   if (!trkbranch)
+       trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
+    trkbranch->SetAddress(&trkl);
+
+    for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
+       if (trackletWord[iTracklet] == 0)
+           continue;
+       trkl->SetTrackletWord(trackletWord[iTracklet]);
+       trkl->SetDetector(30*fSector + 6*fStack + fLayer);
+       trkl->SetROB(fRobPos);
+       trkl->SetMCM(fMcmPos);
+       trkl->SetLabel(trackRef[iTracklet]);
+       trackletTree->Fill();
+    }
+    delete trkl;
+    dl->WriteData("OVERWRITE");
   }
-
-  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;
-  
-*/
+  delete [] adcChannel; 
+  delete [] trackRef;
+  delete digman;
 
   // to be done:
   // error measure for quality of fit (not necessarily needed for the trigger)
   // cluster quality threshold (not yet set)
   // electron probability
-  
+}
+//_____________________________________________________________________________________
+void AliTRDmcmSim::GeneratefZSM1Dim()
+{
+  //
+  // Generate the array fZSM1Dim necessary
+  // for the method ProduceRawStream
+  //
 
-   
+  // Fill the mapping
+  // Supressed zeros indicated by -1 in digits array
+  for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) 
+    {
+      for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
+       {
+         
+         if(fADCF[iadc][it]==-1)  // If is a supressed value
+           {
+             fZSM[iadc][it]=1;
+           }
+         else                    // Not suppressed
+           {
+             fZSM[iadc][it]=0;
+           }
+       }
+    }
+
+  // Make the 1 dim projection
+  for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) 
+    {
+      for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
+       {
+         fZSM1Dim[iadc] &= fZSM[iadc][it];
+       }
+    }
+}
+//_______________________________________________________________________________________
+void AliTRDmcmSim::CopyArrays()
+{
+  //
+  // Initialize filtered data array with raw data
+  // Method added for internal consistency
+  //
+
+  for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) 
+    {
+      for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
+       {
+         fADCF[iadc][it] = fADCR[iadc][it]; 
+       }
+    }
 }
+//_______________________________________________________________________________________
+void AliTRDmcmSim::StartfastZS(Int_t pads, Int_t timebins)
+{
+  //
+  // Initialize just the necessary elements to perform
+  // the zero suppression in the digitizer
+  //
+   
+  fFeeParam  = AliTRDfeeParam::Instance();
+  fSimParam  = AliTRDSimParam::Instance();
+  fNADC      = pads;      
+  fNTimeBin  = timebins; 
 
+  if( fADCR == NULL ) 
+    {
+      fADCR    = new Int_t *[fNADC];
+      fADCF    = new Int_t *[fNADC];
+      fADCT    = new Int_t *[fNADC]; 
+      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]; 
+       fZSM [iadc] = new Int_t[fNTimeBin];
+      }
+    }
 
+  for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) 
+    {
+      for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
+       {
+         fADCR[iadc][it] =  0;
+         fADCF[iadc][it] =  0;
+         fADCT[iadc][it] = -1;  
+         fZSM [iadc][it] =  1;   
+       }
+      fZSM1Dim[iadc] = 1;      
+    }
+  
+  fInitialized = kTRUE;
+}
+//_______________________________________________________________________________________
+void AliTRDmcmSim::FlagDigitsArray(AliTRDarrayADC *tempdigs, Int_t valrow)
+{
+  //
+  // Modify the digits array to flag suppressed values
+  //
 
+  for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) 
+    {
+      for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
+       {
+         if(fZSM[iadc][it]==1)
+           {
+             tempdigs->SetData(valrow,iadc,it,-1);
+           }
+       }
+    }
+}
+//_______________________________________________________________________________________
+void AliTRDmcmSim::RestoreZeros()
+{
+  //
+  // Restore the zero-suppressed values (set as -1) to the value 0
+  //
 
+  for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) 
+    {
+      for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
+       {
+         
+         if(fADCF[iadc][it]==-1)  //if is a supressed zero, reset to zero
+           {
+             fADCF[iadc][it]=0;
+             fADCR[iadc][it]=0;
+           }     
+       }
+    }
 
+}