]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDmcmSim.cxx
AliTRDCalib.cxx .h -> skip all the event if the initialization is not ok for the...
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
index ddb10440f0cb74b7c14cb4e34b3136b016c427cd..c260c53af8f49e81e51e945e366d21962adb124c 100644 (file)
@@ -34,6 +34,7 @@
 #include "TLine.h"
 #include "TRandom.h"
 #include "TClonesArray.h"
+#include "TMath.h"
 
 #include "AliLog.h"
 #include "AliRunLoader.h"
@@ -71,6 +72,7 @@ AliTRDmcmSim::AliTRDmcmSim() :
   fMCMT(NULL),
   fTrackletArray(NULL),
   fZSMap(NULL),
+  fTrklBranchName("mcmtrklbranch"),
   fFeeParam(NULL),
   fTrapConfig(NULL),
   fDigitsManager(NULL),
@@ -86,6 +88,14 @@ AliTRDmcmSim::AliTRDmcmSim() :
   // AliTRDmcmSim default constructor
   // By default, nothing is initialized.
   // It is necessary to issue Init before use.
+
+  for (Int_t iDict = 0; iDict < 3; iDict++)
+    fDict[iDict] = 0x0;
+
+  fFitPtr[0] = 0;
+  fFitPtr[1] = 0;
+  fFitPtr[2] = 0;
+  fFitPtr[3] = 0;
 }
 
 AliTRDmcmSim::~AliTRDmcmSim() 
@@ -467,7 +477,7 @@ void AliTRDmcmSim::Draw(Option_t* const option)
       AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
       Float_t padWidth = 0.635 + 0.03 * (fDetector % 6);
       Float_t offset   = padWidth/256. * ((((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - ((18*4*2 - 18*2 - 3) << 7)); // revert adding offset in FitTracklet
-      Int_t   ndrift   = fTrapConfig->GetDmem(0xc025, fDetector, fRobPos, fMcmPos) >> 5; 
+      Int_t   ndrift   = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrNdrift, fDetector, fRobPos, fMcmPos) >> 5;
       Float_t slope    = trkl->GetdY() * 140e-4 / ndrift; 
 
       Int_t t0 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
@@ -521,7 +531,7 @@ void AliTRDmcmSim::SetData( Int_t adc, Int_t it, Int_t data )
   fADCF[adc][it] = data << fgkAddDigits;
 }
 
-void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
+void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager)
 {
   // Set the ADC data from an AliTRDarrayADC
 
@@ -574,6 +584,64 @@ void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *
   }
 }
 
+void AliTRDmcmSim::SetDataByPad(AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager)
+{
+  // Set the ADC data from an AliTRDarrayADC 
+  // (by pad, to be used during initial reading in simulation)
+
+  if( !CheckInitialized() ) 
+    return;
+
+  fDigitsManager = digitsManager;
+  if (fDigitsManager) {
+    for (Int_t iDict = 0; iDict < 3; iDict++) {
+      AliTRDarrayDictionary *newDict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict);
+      if (fDict[iDict] != 0x0 && newDict != 0x0) {
+        
+        if (fDict[iDict] == newDict)
+          continue;
+
+        fDict[iDict] = newDict;
+        
+        if (fDict[iDict]->GetDim() == 0) {
+          AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
+          continue;
+        }
+        fDict[iDict]->Expand(); 
+      }
+      else {
+        fDict[iDict] = newDict;
+        if (fDict[iDict])
+          fDict[iDict]->Expand();
+      }
+    }
+  }
+
+  if (fNTimeBin != adcArray->GetNtime())
+    SetNTimebins(adcArray->GetNtime());
+  
+  Int_t offset = (fMcmPos % 4 + 1) * 18 + (fRobPos % 2) * 72 + 1;
+
+  for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+    for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
+      Int_t value = -1;
+      Int_t pad = offset - iAdc;
+      if (pad > -1 && pad < 144) 
+       value = adcArray->GetData(GetRow(), offset - iAdc, iTimeBin);
+      //      Int_t value = adcArray->GetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin);
+      if (value < 0 || (offset - iAdc < 1) || (offset - iAdc > 165)) {
+        fADCR[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP) + (fgAddBaseline << fgkAddDigits);
+        fADCF[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
+      }
+      else {
+        fZSMap[iAdc] = 0;
+        fADCR[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits);
+        fADCF[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits);
+      }
+    }
+  }
+}
+
 void AliTRDmcmSim::SetDataPedestal( Int_t adc )
 {
   //
@@ -608,7 +676,7 @@ Bool_t AliTRDmcmSim::GetHit(Int_t index, Int_t &channel, Int_t &timebin, Int_t &
                         (channel << 8) - ypos) 
     * (0.635 + 0.03 * (fDetector % 6))
     / 256.0;
-  label   = fHits[index].fLabel;
+  label   = fHits[index].fLabel[0];
 
   return kTRUE;
 }
@@ -641,6 +709,8 @@ Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t bufSize, UInt_t iEv) co
     return 0;
 
   UInt_t  x;
+  UInt_t  mcmHeader = 0;
+  UInt_t  adcMask = 0;
   Int_t   nw  = 0;  // Number of written words
   Int_t   of  = 0;  // Number of overflowed words
   Int_t   rawVer   = fFeeParam->GetRAWversion();
@@ -655,34 +725,38 @@ Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t bufSize, UInt_t iEv) co
   else 
     adc = fADCF;
   
-  // Produce MCM header
-  x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
-
-  if (nw < bufSize) {
-    buf[nw++] = 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;
+  if( rawVer >= 3 &&
+      (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC15CPUA) & (1 << 13))) { // check for zs flag in TRAP configuration
     for( Int_t iAdc = 0 ; iAdc < fgkNADC ; iAdc++ ) {
       if( ~fZSMap[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
+       adcMask |= (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
 
-    if (nw < bufSize) {
-      buf[nw++] = x;
-    }
-    else {
+    if ((nActiveADC == 0) &&
+       (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC15CPUA) & (1 << 8))) // check for DEH flag in TRAP configuration
+      return 0;
+
+    // assemble adc mask word
+    adcMask |= (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
+  }
+
+  // MCM header
+  mcmHeader = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
+  if (nw < bufSize)
+    buf[nw++] = mcmHeader;
+  else
+    of++;
+
+  // ADC mask
+  if( adcMask != 0 ) {
+    if (nw < bufSize)
+      buf[nw++] = adcMask;
+    else
       of++;
-    }
   }
 
   // Produce ADC data. 3 timebins are packed into one 32 bits word
@@ -1081,7 +1155,7 @@ void AliTRDmcmSim::ZSMapping()
   }
 }
 
-void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label
+void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label[])
 {
   // Add the given hit to the fit register which is lateron used for 
   // the tracklet calculation. 
@@ -1112,7 +1186,9 @@ void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Sh
   fHits[fNHits].fQtot = qtot;
   fHits[fNHits].fYpos = ypos;
   fHits[fNHits].fTimebin = timebin;
-  fHits[fNHits].fLabel = label;
+  fHits[fNHits].fLabel[0] = label[0];
+  fHits[fNHits].fLabel[1] = label[1];
+  fHits[fNHits].fLabel[2] = label[2];
   fNHits++;
 }
 
@@ -1128,7 +1204,7 @@ void AliTRDmcmSim::CalcFitreg()
   
   UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
   Short_t ypos, fromLeft, fromRight, found;
-  UShort_t qTotal[19]; // the last is dummy
+  UShort_t qTotal[19+1]; // the last is dummy
   UShort_t marked[6], qMarked[6], worse1, worse2;
   
   timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS); 
@@ -1282,16 +1358,16 @@ void AliTRDmcmSim::CalcFitreg()
         ypos = 128*(adcLeft - adcRight) / adcCentral;
         if (ypos < 0) ypos = -ypos;
         // make the correction using the position LUT
-        ypos = ypos + fTrapConfig->GetTrapReg((AliTRDtrapConfig::TrapReg_t) (AliTRDtrapConfig::kTPL00 + (ypos & 0x7F))); 
+        ypos = ypos + fTrapConfig->GetTrapReg((AliTRDtrapConfig::TrapReg_t) (AliTRDtrapConfig::kTPL00 + (ypos & 0x7F)),
+                                             fDetector, fRobPos, fMcmPos);
         if (adcLeft > adcRight) ypos = -ypos;
 
-        // label calculation
-        Int_t mcLabel = -1;
+        // label calculation (up to 3)
+        Int_t mcLabel[] = {-1, -1, -1};
         if (fDigitsManager) {
-          Int_t label[9] = { 0 }; // up to 9 different labels possible
-          Int_t count[9] = { 0 };
-          Int_t maxIdx = -1;
-          Int_t maxCount = 0;
+          const Int_t maxLabels = 9;
+          Int_t label[maxLabels] = { 0 }; // up to 9 different labels possible
+          Int_t count[maxLabels] = { 0 };
           Int_t nLabels = 0;
           Int_t padcol[3]; 
           padcol[0] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch);
@@ -1304,26 +1380,29 @@ void AliTRDmcmSim::CalcFitreg()
             for (Int_t iPad = 0; iPad < 3; iPad++) {
               if (padcol[iPad] < 0) 
                 continue;
-              Int_t currLabel = fDict[iDict]->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector);
+              Int_t currLabel = fDict[iDict]->GetData(padrow, padcol[iPad], timebin);
              AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin));
               for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
                 if (currLabel == label[iLabel]) {
                   count[iLabel]++;
-                  if (count[iLabel] > maxCount) {
-                    maxCount = count[iLabel];
-                    maxIdx = iLabel;
-                  }
-                  currLabel = 0;
+                  currLabel = -1;
                   break;
                 }
               } 
-              if (currLabel > 0) {
-                label[nLabels++] = currLabel;
+              if (currLabel >= 0) {
+                label[nLabels] = currLabel;
+               count[nLabels] = 1;
+               nLabels++;
               }
             }
           }
-          if (maxIdx >= 0)
-            mcLabel = label[maxIdx];
+         Int_t index[2*maxLabels];
+         TMath::Sort(maxLabels, count, index);
+         for (Int_t i = 0; i < 3; i++) {
+           if (count[index[i]] <= 0)
+             break;
+           mcLabel[i] = label[index[i]];
+         }
         }
 
         // add the hit to the fitregister
@@ -1444,16 +1523,17 @@ void AliTRDmcmSim::FitTracklet()
   UInt_t scaleY = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 160.0e-4) * shift);
   UInt_t scaleD = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 140.0e-4) * shift);
 
-  Int_t deflCorr = fTrapConfig->GetDmem(0xc022, fDetector, fRobPos, fMcmPos);
-  Int_t ndrift   = fTrapConfig->GetDmem(0xc025, fDetector, fRobPos, fMcmPos); 
+  Int_t deflCorr = (Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCorr, fDetector, fRobPos, fMcmPos);
+  Int_t ndrift   = (Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrNdrift, fDetector, fRobPos, fMcmPos); 
 
   // local variables for calculation
   Long64_t mult, temp, denom; //???
-  UInt_t q0, q1, qTotal;          // charges in the two windows and total charge
+  UInt_t q0, q1, pid;             // charges in the two windows and total charge
   UShort_t nHits;                 // number of hits
   Int_t slope, offset;            // slope and offset of the tracklet
   Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
-  //int32_t SumY2;                // not used in the current TRAP program
+  Int_t sumY2;                // not used in the current TRAP program, now used for error calculation (simulation only)
+  Float_t fitError, fitSlope, fitOffset;
   FitReg_t *fit0, *fit1;          // pointers to relevant fit registers
   
 //  const uint32_t OneDivN[32] = {  // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
@@ -1480,13 +1560,14 @@ void AliTRDmcmSim::FitTracklet()
       nHits   = fit0->fNhits + fit1->fNhits; // number of hits
       sumX    = fit0->fSumX  + fit1->fSumX;
       sumX2   = fit0->fSumX2 + fit1->fSumX2;
-      denom   = nHits*sumX2 - sumX*sumX;
+      denom   = ((Long64_t) nHits)*((Long64_t) sumX2) - ((Long64_t) sumX)*((Long64_t) sumX);
 
       mult    = mult / denom; // exactly like in the TRAP program
       q0      = fit0->fQ0    + fit1->fQ0;
       q1      = fit0->fQ1    + fit1->fQ1;
       sumY    = fit0->fSumY  + fit1->fSumY  + 256*fit1->fNhits;
       sumXY   = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
+      sumY2   = fit0->fSumY2 + fit1->fSumY2 + 512*fit1->fSumY + 256*256*fit1->fNhits;
 
       slope   = nHits*sumXY - sumX * sumY;
       offset  = sumX2*sumY  - sumX * sumXY;
@@ -1515,13 +1596,28 @@ void AliTRDmcmSim::FitTracklet()
 
       AliDebug(5, Form("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i", 
                        fDetector, fRobPos, fMcmPos, slope, 
-                       fTrapConfig->GetDmem(0xc030 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos), 
-                       fTrapConfig->GetDmem(0xc031 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)));
+                       (Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart     + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos), 
+                       (Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)));
+
+      AliDebug(5, Form("Fit sums: x = %i, X = %i, y = %i, Y = %i, Z = %i", 
+                      sumX, sumX2, sumY, sumY2, sumXY));
+
+      fitSlope  = (Float_t) (nHits * sumXY - sumX * sumY) / (nHits * sumX2 - sumX*sumX);
+
+      fitOffset = (Float_t) (sumX2 * sumY - sumX * sumXY) / (nHits * sumX2 - sumX*sumX);
+
+      Float_t sx  = (Float_t) sumX;
+      Float_t sx2 = (Float_t) sumX2;
+      Float_t sy  = (Float_t) sumY;
+      Float_t sy2 = (Float_t) sumY2;
+      Float_t sxy = (Float_t) sumXY;
+      fitError = sy2 - (sx2 * sy*sy - 2 * sx * sxy * sy + nHits * sxy*sxy) / (nHits * sx2 - sx*sx);
+      //fitError = (Float_t) sumY2 - (Float_t) (sumY*sumY) / nHits - fitSlope * ((Float_t) (sumXY - sumX*sumY) / nHits);
 
       Bool_t rejected = kFALSE;
       // deflection range table from DMEM
-      if ((slope < fTrapConfig->GetDmem(0xc030 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)) || 
-          (slope > fTrapConfig->GetDmem(0xc031 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)))
+      if ((slope < ((Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart     + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))) || 
+          (slope > ((Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))))
         rejected = kTRUE;
 
       if (rejected && GetApplyCut())
@@ -1542,25 +1638,23 @@ void AliTRDmcmSim::FitTracklet()
           AliWarning("Overflow in offset");
         offset  = offset & 0x1FFF; // 13 bit
 
-       qTotal = 0; // set to zero as long as no reasonable PID calculation is available 
-                    // before: GetPID(q0/length/fgChargeNorm, q1/length/fgChargeNorm);
+       pid = GetPID(q0 >> fgkAddDigits, q1 >> fgkAddDigits);  // divided by 4 because in simulation there are two additional decimal places
 
-        if (qTotal > 0xff)
-          AliWarning("Overflow in charge");
-        qTotal  = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
+        if (pid > 0xff)
+          AliWarning("Overflow in PID");
+        pid  = pid & 0xFF; // 8 bit, exactly like in the TRAP program
         
         // assemble and store the tracklet word
-        fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
+        fMCMT[cpu] = (pid << 24) | (padrow << 20) | (slope << 13) | offset;
 
         // calculate MC label
-        Int_t mcLabel = -1;
+        Int_t mcLabel[] = { -1, -1, -1};
        Int_t nHits0 = 0;
        Int_t nHits1 = 0;
         if (fDigitsManager) {
-          Int_t label[30] = {0}; // up to 30 different labels possible
-          Int_t count[30] = {0};
-          Int_t maxIdx = -1;
-          Int_t maxCount = 0;
+         const Int_t maxLabels = 30;
+          Int_t label[maxLabels] = {0}; // up to 30 different labels possible
+          Int_t count[maxLabels] = {0};
           Int_t nLabels = 0;
           for (Int_t iHit = 0; iHit < fNHits; iHit++) {
             if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
@@ -1575,24 +1669,29 @@ void AliTRDmcmSim::FitTracklet()
                fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
              nHits1++;
 
-            Int_t currLabel = fHits[iHit].fLabel;
-            for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
-              if (currLabel == label[iLabel]) {
-                count[iLabel]++;
-                if (count[iLabel] > maxCount) {
-                  maxCount = count[iLabel];
-                  maxIdx = iLabel;
-                }
-                currLabel = 0;
-                break;
-              }
-            }
-            if (currLabel > 0) {
-              label[nLabels++] = currLabel;
-            }
-          }
-          if (maxIdx >= 0)
-            mcLabel = label[maxIdx];
+           for (Int_t i = 0; i < 3; i++) {
+             Int_t currLabel = fHits[iHit].fLabel[i];
+             for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
+               if (currLabel == label[iLabel]) {
+                 count[iLabel]++;
+                 currLabel = -1;
+                 break;
+               }
+             }
+             if (currLabel >= 0 && nLabels < maxLabels) {
+               label[nLabels] = currLabel;
+               count[nLabels]++;
+               nLabels++;
+             }
+           }
+         }
+         Int_t index[2*maxLabels];
+         TMath::Sort(maxLabels, count, index);
+         for (Int_t i = 0; i < 3; i++) {
+           if (count[index[i]] <= 0)
+             break;
+           mcLabel[i] = label[index[i]];
+         }
         }
         new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
@@ -1603,6 +1702,36 @@ void AliTRDmcmSim::FitTracklet()
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
+        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetSlope(fitSlope);
+        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetOffset(fitOffset);
+        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetError(TMath::Sqrt(TMath::Abs(fitError)/nHits));
+
+//     // cluster information
+//     Float_t *res = new Float_t[nHits];
+//     Float_t *qtot = new Float_t[nHits];
+//     Int_t nCls = 0;
+//     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
+//       // check if hit contributes
+//       if (fHits[iHit].fChannel == fFitPtr[cpu]) {
+//         res[nCls] = fHits[iHit].fYpos - (fitSlope * fHits[iHit].fTimebin + fitOffset);
+//         qtot[nCls] = fHits[iHit].fQtot;
+//         nCls++;
+//       }
+//       else if (fHits[iHit].fChannel == fFitPtr[cpu] + 1) {
+//         res[nCls] = fHits[iHit].fYpos + 256 - (fitSlope * fHits[iHit].fTimebin + fitOffset);
+//         qtot[nCls] = fHits[iHit].fQtot;
+//         nCls++;
+//       }
+//     }
+//        ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetClusters(res, qtot, nCls);
+//     delete [] res;
+//     delete [] qtot;
+
+       if (fitError < 0)
+         AliError(Form("Strange fit error: %f from Sx: %i, Sy: %i, Sxy: %i, Sx2: %i, Sy2: %i, nHits: %i",
+                       fitError, sumX, sumY, sumXY, sumX2, sumY2, nHits));
+       AliDebug(3, Form("fit slope: %f, offset: %f, error: %f", 
+                        fitSlope, fitOffset, TMath::Sqrt(TMath::Abs(fitError)/nHits)));
       }
     }
   }
@@ -1651,9 +1780,9 @@ Bool_t AliTRDmcmSim::StoreTracklets()
   }
   
   AliTRDtrackletMCM *trkl = 0x0;
-  TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
+  TBranch *trkbranch = trackletTree->GetBranch(fTrklBranchName.Data());
   if (!trkbranch)
-    trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
+    trkbranch = trackletTree->Branch(fTrklBranchName.Data(), "AliTRDtrackletMCM", &trkl, 32000);
   
   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
     trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
@@ -1683,6 +1812,11 @@ void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
           digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, -1);
         }
       }
+      else if (iAdc < 2 || iAdc == 20) {
+        for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+          digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, (fADCR[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
+        }
+      }
     }
   }
   else {
@@ -1701,6 +1835,81 @@ void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
   }
 }
 
+
+// ******************************
+// PID section
+//
+// Memory area for the LUT: 0xC100 to 0xC3FF
+//
+// The addresses for the parameters (the order is optimized for maximum calculation speed in the MCMs):
+// 0xC028: cor1
+// 0xC029: nBins(sF)
+// 0xC02A: cor0
+// 0xC02B: TableLength
+// Defined in AliTRDtrapConfig.h
+//
+// The algorithm implemented in the TRAP program of the MCMs (Venelin Angelov)
+//  1) set the read pointer to the beginning of the Parameters in DMEM
+//  2) shift right the FitReg with the Q0 + (Q1 << 16) to get Q1
+//  3) read cor1 with rpointer++
+//  4) start cor1*Q1
+//  5) read nBins with rpointer++
+//  6) start nBins*cor1*Q1
+//  7) read cor0 with rpointer++
+//  8) swap hi-low parts in FitReg, now is Q1 + (Q0 << 16)
+//  9) shift right to get Q0
+// 10) start cor0*Q0
+// 11) read TableLength
+// 12) compare cor0*Q0 with nBins
+// 13) if >=, clip cor0*Q0 to nBins-1
+// 14) add cor0*Q0 to nBins*cor1*Q1
+// 15) compare the result with TableLength
+// 16) if >=, clip to TableLength-1
+// 17) read from the LUT 8 bits
+
+
+Int_t AliTRDmcmSim::GetPID(Int_t q0, Int_t q1)
+{
+  // return PID calculated from charges accumulated in two time windows
+
+   ULong64_t addrQ0;
+   ULong64_t addr;
+
+   UInt_t nBinsQ0 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTnbins);  // number of bins in q0 / 4 !!
+   UInt_t pidTotalSize = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength);
+   if(nBinsQ0==0 || pidTotalSize==0)  // make sure we don't run into trouble if one of the values is not configured
+      return 0;
+
+   ULong_t corrQ0 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTcor0, fDetector, fRobPos, fMcmPos);
+   ULong_t corrQ1 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTcor1, fDetector, fRobPos, fMcmPos);
+   if(corrQ0==0 || corrQ1==0)  // make sure we don't run into trouble if one of the values is not configured
+      return 0;
+
+   addrQ0 = corrQ0;
+   addrQ0 = (((addrQ0*q0)>>16)>>16); // because addrQ0 = (q0 * corrQ0) >> 32; does not work for unknown reasons
+
+   if(addrQ0 >= nBinsQ0) {  // check for overflow
+      AliDebug(5,Form("Overflow in q0: %llu/4 is bigger then %u", addrQ0, nBinsQ0));
+      addrQ0 = nBinsQ0 -1;
+   } 
+
+   addr = corrQ1;
+   addr = (((addr*q1)>>16)>>16);
+   addr = addrQ0 + nBinsQ0*addr; // because addr = addrQ0 + nBinsQ0* (((corrQ1*q1)>>32); does not work
+
+   if(addr >= pidTotalSize) {
+      AliDebug(5,Form("Overflow in q1. Address %llu/4 is bigger then %u", addr, pidTotalSize));
+      addr = pidTotalSize -1;
+   } 
+
+   // For a LUT with 11 input and 8 output bits, the first memory address is set to  LUT[0] | (LUT[1] << 8) | (LUT[2] << 16) | (LUT[3] << 24)
+   // and so on
+   UInt_t result = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTStart+(addr/4));
+   return (result>>((addr%4)*8)) & 0xFF;
+}
+
+
+
 // help functions, to be cleaned up
 
 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
@@ -1727,8 +1936,8 @@ UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
 
 void AliTRDmcmSim::Sort2(UShort_t  idx1i, UShort_t  idx2i, \
                             UShort_t  val1i, UShort_t  val2i, \
-                            UShort_t *idx1o, UShort_t *idx2o, \
-                            UShort_t *val1o, UShort_t *val2o) const
+                            UShort_t * const idx1o, UShort_t * const idx2o, \
+                            UShort_t * const val1o, UShort_t * const val2o) const
 {
   // sorting for tracklet selection
 
@@ -1750,8 +1959,8 @@ void AliTRDmcmSim::Sort2(UShort_t  idx1i, UShort_t  idx2i, \
 
 void AliTRDmcmSim::Sort3(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, \
                             UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, \
-                            UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \
-                            UShort_t *val1o, UShort_t *val2o, UShort_t *val3o)
+                            UShort_t * const idx1o, UShort_t * const idx2o, UShort_t * const idx3o, \
+                            UShort_t * const val1o, UShort_t * const val2o, UShort_t * const val3o)
 {
   // sorting for tracklet selection
 
@@ -1826,8 +2035,8 @@ void AliTRDmcmSim::Sort3(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, \
 
 void AliTRDmcmSim::Sort6To4(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
                                UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
-                               UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \
-                               UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o)
+                               UShort_t * const idx1o, UShort_t * const idx2o, UShort_t * const idx3o, UShort_t * const idx4o, \
+                               UShort_t * const val1o, UShort_t * const val2o, UShort_t * const val3o, UShort_t * const val4o)
 {
   // sorting for tracklet selection
 
@@ -1854,7 +2063,7 @@ void AliTRDmcmSim::Sort6To4(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, U
 
 void AliTRDmcmSim::Sort6To2Worst(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
                                     UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
-                                    UShort_t *idx5o, UShort_t *idx6o)
+                                    UShort_t * const idx5o, UShort_t * const idx6o)
 {
   // sorting for tracklet selection
 
@@ -1985,3 +2194,201 @@ ostream& operator<<(ostream& os, const AliTRDmcmSim& mcm)
 
   return os;
 }
+
+
+void AliTRDmcmSim::PrintFitRegXml(ostream& os) const
+{
+  // print fit registres in XML format
+
+   bool tracklet=false;
+
+  for (Int_t cpu = 0; cpu < 4; cpu++) {
+     if(fFitPtr[cpu] != 31)
+       tracklet=true;
+  }
+
+  if(tracklet==true) {
+     os << "<nginject>" << std::endl;
+     os << "<ack roc=\""<< fDetector <<  "\" cmndid=\"0\">" << std::endl;
+     os << "<dmem-readout>" << std::endl;
+     os << "<d det=\"" << fDetector << "\">" << std::endl;
+     os << " <ro-board rob=\"" << fRobPos << "\">" << std::endl;
+     os << "  <m mcm=\"" << fMcmPos << "\">" << std::endl;
+     
+     for(int cpu=0; cpu<4; cpu++) {
+       os << "   <c cpu=\"" << cpu << "\">" << std::endl;
+       if(fFitPtr[cpu] != 31) {
+          for(int adcch=fFitPtr[cpu]; adcch<fFitPtr[cpu]+2; adcch++) {
+             os << "    <ch chnr=\"" << adcch << "\">"<< std::endl;
+             os << "     <hits>"   << fFitReg[adcch].fNhits << "</hits>"<< std::endl;
+             os << "     <q0>"     << fFitReg[adcch].fQ0/4 << "</q0>"<< std::endl;    // divided by 4 because in simulation we have 2 additional decimal places
+             os << "     <q1>"     << fFitReg[adcch].fQ1/4 << "</q1>"<< std::endl;    // in the output 
+             os << "     <sumx>"   << fFitReg[adcch].fSumX << "</sumx>"<< std::endl;
+             os << "     <sumxsq>" << fFitReg[adcch].fSumX2 << "</sumxsq>"<< std::endl;
+             os << "     <sumy>"   << fFitReg[adcch].fSumY << "</sumy>"<< std::endl;
+             os << "     <sumysq>" << fFitReg[adcch].fSumY2 << "</sumysq>"<< std::endl;
+             os << "     <sumxy>"  << fFitReg[adcch].fSumXY << "</sumxy>"<< std::endl;
+             os << "    </ch>" << std::endl;
+          }
+       }
+       os << "      </c>" << std::endl;
+     }
+     os << "    </m>" << std::endl;
+     os << "  </ro-board>" << std::endl;
+     os << "</d>" << std::endl;
+     os << "</dmem-readout>" << std::endl;
+     os << "</ack>" << std::endl;
+     os << "</nginject>" << std::endl;
+  }
+}
+
+
+void AliTRDmcmSim::PrintTrackletsXml(ostream& os) const
+{
+  // print tracklets in XML format
+
+   os << "<nginject>" << std::endl;
+   os << "<ack roc=\""<< fDetector <<  "\" cmndid=\"0\">" << std::endl;
+   os << "<dmem-readout>" << std::endl;
+   os << "<d det=\"" << fDetector << "\">" << std::endl;
+   os << "  <ro-board rob=\"" << fRobPos << "\">" << std::endl;
+   os << "    <m mcm=\"" << fMcmPos << "\">" << std::endl;
+
+   Int_t pid, padrow, slope, offset;
+   for(Int_t cpu=0; cpu<4; cpu++) {
+      if(fMCMT[cpu] == 0x10001000) {
+        pid=-1;
+        padrow=-1;
+        slope=-1;
+        offset=-1;
+      }
+      else {
+        pid    = (fMCMT[cpu] & 0xFF000000) >> 24;
+        padrow = (fMCMT[cpu] & 0xF00000  ) >> 20;
+        slope  = (fMCMT[cpu] & 0xFE000   ) >> 13;
+        offset = (fMCMT[cpu] & 0x1FFF    ) ;
+
+      }
+      os << "      <trk> <pid>" << pid << "</pid>" << " <padrow>" << padrow << "</padrow>" 
+        << " <slope>" << slope << "</slope>" << " <offset>" << offset << "</offset>" << "</trk>" << std::endl;
+   }
+
+   os << "    </m>" << std::endl;
+   os << "  </ro-board>" << std::endl;
+   os << "</d>" << std::endl;
+   os << "</dmem-readout>" << std::endl;
+   os << "</ack>" << std::endl;
+   os << "</nginject>" << std::endl;
+}
+
+
+void AliTRDmcmSim::PrintAdcDatHuman(ostream& os) const
+{
+  // print ADC data in human-readable format
+
+   os << "MCM " << fMcmPos << " on ROB " << fRobPos << 
+      " in detector " << fDetector << std::endl;
+    
+   os << "----- Unfiltered ADC data (10 bit) -----" << std::endl;
+   os << "ch    ";
+   for (Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) 
+      os << std::setw(5) << iChannel;
+   os << std::endl;
+   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+      os << "tb " << std::setw(2) << iTimeBin << ":";
+      for (Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) {
+        os << std::setw(5) << (fADCR[iChannel][iTimeBin] >> fgkAddDigits);
+      }
+      os << std::endl;
+   }
+    
+   os << "----- Filtered ADC data (10+2 bit) -----" << std::endl;
+   os << "ch    ";
+   for (Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) 
+      os << std::setw(4) << iChannel
+         << ((~fZSMap[iChannel] != 0) ? "!" : " ");
+   os << std::endl;
+   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+      os << "tb " << std::setw(2) << iTimeBin << ":";
+      for (Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) {
+        os << std::setw(4) << (fADCF[iChannel][iTimeBin])
+           << (((fZSMap[iChannel] & (1 << iTimeBin)) == 0) ? "!" : " ");
+      }
+      os << std::endl;
+   }
+}
+
+
+void AliTRDmcmSim::PrintAdcDatXml(ostream& os) const
+{
+  // print ADC data in XML format 
+
+   os << "<nginject>" << std::endl;
+   os << "<ack roc=\""<< fDetector <<  "\" cmndid=\"0\">" << std::endl;
+   os << "<dmem-readout>" << std::endl;
+   os << "<d det=\"" << fDetector << "\">" << std::endl;
+   os << " <ro-board rob=\"" << fRobPos << "\">" << std::endl;
+   os << "  <m mcm=\"" << fMcmPos << "\">" << std::endl;
+
+    for(Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) {
+       os << "   <ch chnr=\"" << iChannel << "\">" << std::endl;
+       for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+         os << "<tb>" << fADCF[iChannel][iTimeBin]/4 << "</tb>";
+       }
+       os << "   </ch>" << std::endl;
+    }
+
+   os << "  </m>" << std::endl;
+   os << " </ro-board>" << std::endl;
+   os << "</d>" << std::endl;
+   os << "</dmem-readout>" << std::endl;
+   os << "</ack>" << std::endl;
+   os << "</nginject>" << std::endl;
+}
+
+
+
+void AliTRDmcmSim::PrintAdcDatDatx(ostream& os, Bool_t broadcast) const
+{
+  // print ADC data in datx format (to send to FEE)
+
+   fTrapConfig->PrintDatx(os, 2602, 1, 0, 127);  // command to enable the ADC clock - necessary to write ADC values to MCM
+   os << std::endl;
+
+   Int_t addrOffset = 0x2000;
+   Int_t addrStep   = 0x80;
+   Int_t addrOffsetEBSIA = 0x20;
+    
+   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+      for (Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) {
+        if(broadcast==kFALSE)
+           fTrapConfig->PrintDatx(os, addrOffset+iChannel*addrStep+addrOffsetEBSIA+iTimeBin, (fADCF[iChannel][iTimeBin]/4), GetRobPos(),  GetMcmPos());
+        else
+           fTrapConfig->PrintDatx(os, addrOffset+iChannel*addrStep+addrOffsetEBSIA+iTimeBin, (fADCF[iChannel][iTimeBin]/4), 0, 127);
+      }
+      os << std::endl;
+   }
+}
+
+
+void AliTRDmcmSim::PrintPidLutHuman()
+{
+  // print PID LUT in human readable format
+
+   UInt_t result;
+
+   UInt_t addrEnd = AliTRDtrapConfig::fgkDmemAddrLUTStart + fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength)/4; // /4 because each addr contains 4 values
+   UInt_t nBinsQ0 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTnbins);
+
+   std::cout << "nBinsQ0: " << nBinsQ0 << std::endl;
+   std::cout << "LUT table length: " << fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength) << std::endl;
+   for(UInt_t addr=AliTRDtrapConfig::fgkDmemAddrLUTStart; addr< addrEnd; addr++) {
+      result = fTrapConfig->GetDmemUnsigned(addr);
+      std::cout << addr << " # x: " << ((addr-AliTRDtrapConfig::fgkDmemAddrLUTStart)%((nBinsQ0)/4))*4 << ", y: " <<(addr-AliTRDtrapConfig::fgkDmemAddrLUTStart)/(nBinsQ0/4)
+               << "  #  " <<((result>>0)&0xFF)
+               << " | "  << ((result>>8)&0xFF)
+               << " | "  << ((result>>16)&0xFF)
+               << " | "  << ((result>>24)&0xFF) << std::endl;
+   }
+}