TRD/AliTRDmcmSim.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  TRD MCM (Multi Chip Module) simulator                                    //
  21 //  which simulated the TRAP processing after the AD-conversion              //
  22 //  The relevant parameters (i.e. configuration registers of the TRAP        //
  23 //  configuration are taken from AliTRDtrapConfig.                           //
  24 //                                                                           //
  25 ///////////////////////////////////////////////////////////////////////////////
  26
  27 #include <fstream>  // needed for raw data dump
  28
  29 #include <TCanvas.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TGraph.h>
  33 #include <TLine.h>
  34 #include <TMath.h>
  35 #include <TRandom.h>
  36 #include <TClonesArray.h>
  37
  38 #include "AliLog.h"
  39 #include "AliRun.h"
  40 #include "AliRunLoader.h"
  41 #include "AliLoader.h"
  42 #include "AliTRDdigit.h"
  43
  44 #include "AliTRDfeeParam.h"
  45 #include "AliTRDtrapConfig.h"
  46 #include "AliTRDSimParam.h"
  47 #include "AliTRDgeometry.h"
  48 #include "AliTRDcalibDB.h"
  49 #include "AliTRDdigitsManager.h"
  50 #include "AliTRDarrayADC.h"
  51 #include "AliTRDarrayDictionary.h"
  52 #include "AliTRDpadPlane.h"
  53 #include "AliTRDtrackletMCM.h"
  54 #include "AliTRDmcmSim.h"
  55
  56 #include "AliMagF.h"
  57 #include "TGeoGlobalMagField.h"
  58
  59 ClassImp(AliTRDmcmSim)
  60
  61 Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
  62
  63 Float_t AliTRDmcmSim::fgChargeNorm = 65000.;
  64 Int_t AliTRDmcmSim::fgAddBaseline = 10;
  65
  66 Int_t AliTRDmcmSim::fgPidNBinsQ0 = 40;
  67 Int_t AliTRDmcmSim::fgPidNBinsQ1 = 50;
  68 Bool_t AliTRDmcmSim::fgPidLutDelete = kFALSE;
  69 Int_t AliTRDmcmSim::fgPidLutDefault[40][50] = {
  70   { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0  },
  71   { 0, 9, 6, 12, 29, 53, 76, 94, 107, 116, 122, 126, 128, 129, 129, 129, 128, 127, 126, 124, 122, 120, 117, 115, 112, 109, 107, 104, 101, 99, 96, 94, 91, 89, 87, 85, 83, 81, 79, 78, 77, 75, 74, 73, 72, 72, 71, 71, 70, 70  },
  72   { 0, 14, 8, 17, 37, 66, 94, 116, 131, 140, 146, 150, 152, 153, 153, 152, 150, 148, 145, 143, 139, 136, 132, 129, 125, 121, 118, 114, 110, 107, 104, 101, 98, 95, 93, 91, 89, 87, 85, 83, 82, 81, 80, 79, 78, 77, 77, 76, 76, 75  },
  73   { 0, 33, 19, 34, 69, 112, 145, 167, 181, 189, 194, 196, 197, 197, 196, 194, 191, 188, 184, 180, 175, 170, 165, 159, 154, 148, 143, 137, 132, 127, 123, 118, 114, 111, 107, 104, 101, 99, 96, 94, 92, 91, 89, 88, 87, 86, 85, 85, 84, 84  },
  74   { 0, 82, 52, 83, 136, 180, 205, 218, 226, 230, 232, 233, 233, 233, 232, 230, 228, 226, 223, 219, 215, 210, 205, 199, 193, 187, 180, 173, 167, 160, 154, 148, 142, 136, 131, 127, 122, 119, 115, 112, 109, 106, 104, 102, 100, 99, 97, 96, 95, 94  },
  75   { 0, 132, 96, 136, 185, 216, 231, 238, 242, 244, 245, 245, 245, 245, 245, 244, 243, 242, 240, 238, 236, 233, 230, 226, 222, 217, 212, 206, 200, 193, 187, 180, 173, 167, 161, 155, 149, 144, 139, 134, 130, 126, 123, 120, 117, 114, 112, 110, 108, 107  },
  76   { 0, 153, 120, 160, 203, 227, 238, 243, 246, 247, 248, 249, 249, 249, 248, 248, 247, 246, 245, 244, 243, 241, 239, 237, 234, 231, 228, 224, 219, 215, 209, 204, 198, 192, 186, 180, 174, 168, 163, 157, 152, 147, 143, 139, 135, 131, 128, 125, 123, 120  },
  77   { 0, 156, 128, 166, 207, 229, 239, 244, 247, 248, 249, 249, 249, 249, 249, 249, 248, 247, 247, 246, 244, 243, 242, 240, 238, 236, 233, 230, 227, 224, 220, 216, 212, 207, 202, 197, 192, 187, 181, 176, 171, 166, 161, 156, 152, 148, 144, 140, 137, 134  },
  78   { 0, 152, 128, 166, 206, 228, 239, 244, 246, 248, 249, 249, 249, 249, 249, 248, 248, 247, 246, 245, 244, 243, 241, 240, 238, 236, 234, 232, 229, 226, 224, 220, 217, 214, 210, 206, 202, 197, 193, 188, 184, 179, 174, 170, 166, 161, 157, 153, 150, 146  },
  79   { 0, 146, 126, 164, 203, 226, 237, 243, 246, 247, 248, 248, 248, 248, 248, 247, 247, 246, 245, 244, 242, 241, 239, 238, 236, 234, 232, 230, 227, 225, 223, 220, 217, 215, 212, 209, 205, 202, 199, 195, 191, 187, 183, 179, 175, 171, 168, 164, 160, 156  },
  80   { 0, 140, 123, 160, 200, 224, 235, 241, 244, 246, 247, 247, 247, 247, 247, 246, 245, 244, 243, 242, 240, 238, 237, 235, 233, 230, 228, 226, 224, 221, 219, 217, 215, 212, 210, 207, 205, 202, 200, 197, 194, 191, 188, 184, 181, 178, 174, 171, 168, 164  },
  81   { 0, 133, 119, 156, 196, 220, 233, 239, 243, 245, 245, 246, 246, 246, 245, 244, 243, 242, 241, 239, 237, 235, 233, 231, 229, 226, 224, 221, 219, 216, 214, 212, 210, 208, 206, 204, 202, 199, 197, 195, 193, 191, 188, 186, 183, 181, 178, 175, 172, 169  },
  82   { 0, 127, 115, 152, 192, 217, 230, 237, 241, 243, 244, 244, 244, 244, 243, 242, 241, 240, 238, 236, 234, 232, 229, 227, 224, 221, 218, 216, 213, 210, 208, 206, 203, 201, 200, 198, 196, 194, 193, 191, 190, 188, 186, 185, 183, 181, 179, 177, 174, 172  },
  83   { 0, 121, 111, 147, 187, 213, 227, 235, 239, 241, 242, 243, 243, 242, 241, 240, 239, 237, 236, 233, 231, 228, 225, 222, 219, 216, 213, 210, 207, 204, 201, 199, 196, 194, 192, 191, 189, 188, 187, 185, 184, 183, 182, 181, 180, 178, 177, 176, 174, 172  },
  84   { 0, 116, 107, 142, 181, 209, 224, 232, 237, 239, 240, 241, 241, 240, 239, 238, 237, 235, 233, 230, 227, 224, 221, 218, 214, 211, 207, 204, 200, 197, 194, 191, 189, 187, 185, 183, 182, 180, 179, 178, 178, 177, 176, 175, 175, 174, 173, 172, 172, 170  },
  85   { 0, 112, 103, 136, 176, 204, 220, 229, 234, 237, 238, 239, 239, 238, 237, 236, 234, 232, 230, 227, 224, 221, 217, 213, 209, 205, 201, 198, 194, 190, 187, 184, 181, 179, 177, 175, 174, 172, 171, 171, 170, 169, 169, 169, 168, 168, 168, 168, 167, 167  },
  86   { 0, 107, 99, 131, 170, 199, 216, 226, 231, 234, 236, 237, 237, 236, 235, 234, 232, 230, 227, 224, 221, 217, 213, 209, 205, 200, 196, 192, 188, 184, 180, 177, 174, 172, 169, 167, 166, 164, 163, 162, 162, 161, 161, 161, 161, 161, 161, 162, 162, 162  },
  87   { 0, 104, 94, 125, 164, 193, 212, 222, 228, 232, 233, 234, 234, 234, 233, 231, 229, 227, 224, 221, 218, 214, 210, 205, 201, 196, 191, 187, 182, 178, 174, 171, 168, 165, 162, 160, 158, 157, 155, 154, 154, 153, 153, 153, 153, 154, 154, 154, 155, 155  },
  88   { 0, 100, 90, 119, 157, 188, 207, 219, 225, 229, 231, 232, 232, 231, 230, 229, 227, 224, 222, 218, 215, 211, 206, 202, 197, 192, 187, 182, 178, 173, 169, 165, 162, 158, 156, 153, 151, 149, 148, 147, 146, 146, 145, 145, 145, 146, 146, 147, 148, 148  },
  89   { 0, 97, 86, 113, 150, 182, 202, 215, 222, 226, 228, 229, 229, 229, 228, 226, 224, 222, 219, 216, 212, 208, 203, 199, 194, 188, 183, 178, 173, 169, 164, 160, 156, 153, 150, 147, 145, 143, 141, 140, 139, 138, 138, 138, 138, 138, 139, 139, 140, 141  },
  90   { 0, 94, 82, 107, 144, 176, 197, 210, 218, 223, 225, 227, 227, 227, 226, 224, 222, 220, 217, 213, 209, 205, 201, 196, 191, 186, 180, 175, 170, 165, 160, 156, 152, 148, 145, 142, 139, 137, 135, 134, 132, 131, 131, 131, 131, 131, 131, 132, 132, 133  },
  91   { 0, 92, 78, 102, 137, 169, 192, 206, 215, 220, 223, 224, 224, 224, 223, 222, 220, 217, 215, 211, 207, 203, 199, 194, 188, 183, 178, 172, 167, 162, 157, 152, 148, 144, 140, 137, 134, 132, 130, 128, 127, 125, 125, 124, 124, 124, 124, 125, 125, 126  },
  92   { 0, 90, 75, 96, 131, 163, 187, 202, 211, 216, 220, 221, 222, 222, 221, 220, 218, 215, 212, 209, 205, 201, 197, 192, 187, 181, 176, 170, 165, 159, 154, 149, 145, 141, 137, 133, 130, 128, 125, 123, 122, 120, 119, 118, 118, 118, 118, 118, 119, 119  },
  93   { 0, 88, 71, 91, 124, 157, 181, 197, 207, 213, 217, 219, 219, 219, 219, 217, 216, 213, 211, 207, 204, 200, 195, 190, 185, 180, 174, 169, 163, 158, 152, 147, 142, 138, 134, 130, 127, 124, 121, 119, 117, 116, 114, 114, 113, 112, 112, 112, 112, 113  },
  94   { 0, 87, 68, 86, 118, 151, 176, 192, 203, 210, 214, 216, 217, 217, 217, 215, 214, 212, 209, 206, 202, 198, 194, 189, 184, 179, 173, 167, 162, 156, 151, 146, 141, 136, 132, 128, 124, 121, 118, 116, 114, 112, 110, 109, 108, 108, 107, 107, 107, 107  },
  95   { 0, 85, 65, 81, 112, 144, 170, 188, 199, 206, 211, 213, 214, 215, 214, 213, 212, 210, 207, 204, 201, 197, 193, 188, 183, 178, 172, 167, 161, 155, 150, 145, 140, 135, 130, 126, 122, 119, 116, 113, 111, 109, 107, 106, 105, 104, 103, 103, 102, 102  },
  96   { 0, 84, 62, 77, 106, 138, 165, 183, 195, 203, 208, 210, 212, 212, 212, 211, 210, 208, 206, 203, 200, 196, 192, 187, 183, 177, 172, 166, 161, 155, 150, 144, 139, 134, 129, 125, 121, 117, 114, 111, 109, 106, 104, 103, 101, 100, 99, 99, 98, 98  },
  97   { 0, 84, 60, 73, 101, 133, 159, 178, 191, 199, 204, 208, 209, 210, 210, 209, 208, 206, 204, 202, 199, 195, 191, 187, 182, 177, 172, 166, 161, 155, 150, 144, 139, 134, 129, 124, 120, 116, 113, 110, 107, 104, 102, 100, 99, 98, 96, 96, 95, 95  },
  98   { 0, 83, 58, 69, 96, 127, 154, 174, 187, 196, 201, 205, 207, 208, 208, 207, 206, 205, 203, 200, 197, 194, 190, 186, 182, 177, 172, 167, 161, 156, 150, 145, 139, 134, 129, 124, 120, 116, 112, 109, 106, 103, 101, 99, 97, 95, 94, 93, 92, 92  },
  99   { 0, 82, 56, 66, 91, 121, 149, 169, 183, 192, 198, 202, 204, 206, 206, 206, 205, 203, 201, 199, 196, 193, 190, 186, 182, 177, 172, 167, 162, 156, 151, 145, 140, 135, 129, 125, 120, 116, 112, 108, 105, 102, 100, 97, 95, 94, 92, 91, 90, 89  },
 100   { 0, 82, 54, 62, 86, 116, 144, 165, 179, 189, 195, 199, 202, 203, 204, 204, 203, 202, 200, 198, 196, 193, 189, 186, 182, 177, 173, 168, 163, 157, 152, 146, 141, 136, 130, 125, 121, 116, 112, 108, 105, 102, 99, 96, 94, 92, 91, 89, 88, 87  },
 101   { 0, 82, 52, 59, 82, 111, 139, 160, 175, 185, 192, 197, 200, 201, 202, 202, 201, 200, 199, 197, 195, 192, 189, 186, 182, 178, 173, 168, 163, 158, 153, 148, 142, 137, 132, 127, 122, 117, 113, 109, 105, 102, 99, 96, 94, 92, 90, 88, 87, 85  },
 102   { 0, 82, 50, 56, 78, 106, 134, 156, 171, 182, 189, 194, 197, 199, 200, 200, 200, 199, 198, 196, 194, 191, 188, 185, 182, 178, 174, 169, 164, 159, 154, 149, 144, 138, 133, 128, 123, 118, 114, 110, 106, 102, 99, 96, 93, 91, 89, 87, 86, 84  },
 103   { 0, 82, 49, 54, 74, 102, 129, 151, 167, 179, 186, 191, 195, 197, 198, 198, 198, 197, 196, 195, 193, 191, 188, 185, 182, 178, 174, 170, 165, 161, 156, 151, 145, 140, 135, 130, 125, 120, 115, 111, 107, 103, 100, 97, 94, 91, 89, 87, 85, 83  },
 104   { 0, 82, 47, 51, 70, 97, 124, 147, 164, 175, 183, 189, 192, 195, 196, 197, 197, 196, 195, 194, 192, 190, 188, 185, 182, 178, 175, 171, 166, 162, 157, 152, 147, 142, 137, 132, 127, 122, 117, 112, 108, 104, 101, 97, 94, 91, 89, 87, 85, 83  },
 105   { 0, 83, 46, 49, 67, 93, 120, 143, 160, 172, 180, 186, 190, 192, 194, 195, 195, 195, 194, 193, 191, 189, 187, 185, 182, 179, 175, 172, 167, 163, 159, 154, 149, 144, 139, 134, 129, 124, 119, 114, 110, 106, 102, 98, 95, 92, 89, 87, 85, 83  },
 106   { 0, 83, 45, 47, 64, 89, 116, 139, 156, 169, 177, 184, 188, 190, 192, 193, 193, 193, 193, 192, 190, 189, 187, 184, 182, 179, 176, 172, 168, 164, 160, 156, 151, 146, 141, 136, 131, 126, 121, 116, 112, 108, 104, 100, 96, 93, 90, 88, 85, 83  },
 107   { 0, 84, 44, 45, 61, 85, 111, 134, 152, 165, 175, 181, 185, 188, 190, 191, 192, 192, 191, 191, 189, 188, 186, 184, 182, 179, 176, 173, 169, 166, 162, 157, 153, 148, 143, 138, 133, 128, 124, 119, 114, 110, 106, 102, 98, 95, 91, 89, 86, 84  },
 108   { 0, 85, 43, 43, 58, 81, 107, 131, 149, 162, 172, 178, 183, 186, 188, 190, 190, 190, 190, 189, 188, 187, 186, 184, 182, 179, 176, 173, 170, 167, 163, 159, 155, 150, 145, 141, 136, 131, 126, 121, 117, 112, 108, 104, 100, 96, 93, 90, 87, 85  },
 109   { 0, 85, 42, 41, 55, 78, 103, 127, 145, 159, 169, 176, 181, 184, 186, 188, 189, 189, 189, 188, 188, 186, 185, 183, 181, 179, 177, 174, 171, 168, 164, 160, 156, 152, 148, 143, 138, 134, 129, 124, 119, 115, 110, 106, 102, 98, 95, 91, 88, 86  }
 110 };
 111
 112 Int_t (*AliTRDmcmSim::fgPidLut) = *fgPidLutDefault;
 113
 114 //_____________________________________________________________________________
 115 AliTRDmcmSim::AliTRDmcmSim() : TObject()
 116   ,fInitialized(kFALSE)
 117   ,fMaxTracklets(-1)
 118   ,fDetector(-1)
 119   ,fRobPos(-1)
 120   ,fMcmPos(-1)
 121   ,fRow (-1)
 122   ,fNADC(-1)
 123   ,fNTimeBin(-1)
 124   ,fADCR(NULL)
 125   ,fADCF(NULL)
 126   ,fMCMT(NULL)
 127   ,fTrackletArray(NULL)
 128   ,fZSM(NULL)
 129   ,fZSM1Dim(NULL)
 130   ,fFeeParam(NULL)
 131   ,fTrapConfig(NULL)
 132   ,fSimParam(NULL)
 133   ,fCommonParam(NULL)
 134   ,fCal(NULL)
 135   ,fGeo(NULL)
 136   ,fDigitsManager(NULL)
 137   ,fPedAcc(NULL)
 138   ,fGainCounterA(NULL)
 139   ,fGainCounterB(NULL)
 140   ,fTailAmplLong(NULL)
 141   ,fTailAmplShort(NULL)
 142   ,fNHits(0)
 143   ,fFitReg(NULL)
 144 {
 145   //
 146   // AliTRDmcmSim default constructor
 147   // By default, nothing is initialized.
 148   // It is necessary to issue Init before use.
 149 }
 150
 151 AliTRDmcmSim::~AliTRDmcmSim()
 152 {
 153   //
 154   // AliTRDmcmSim destructor
 155   //
 156
 157   if(fInitialized) {
 158     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 159       delete [] fADCR[iadc];
 160       delete [] fADCF[iadc];
 161       delete [] fZSM [iadc];
 162     }
 163     delete [] fADCR;
 164     delete [] fADCF;
 165     delete [] fZSM;
 166     delete [] fZSM1Dim;
 167     delete [] fMCMT;
 168
 169     delete [] fPedAcc;
 170     delete [] fGainCounterA;
 171     delete [] fGainCounterB;
 172     delete [] fTailAmplLong;
 173     delete [] fTailAmplShort;
 174     delete [] fFitReg;
 175
 176     fTrackletArray->Delete();
 177     delete fTrackletArray;
 178     delete fGeo;
 179   }
 180 }
 181
 182 void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
 183 {
 184   //
 185   // Initialize the class with new geometry information
 186   // fADC array will be reused with filled by zero
 187   //
 188
 189   if (!fInitialized) {
 190     fFeeParam      = AliTRDfeeParam::Instance();
 191     fTrapConfig    = AliTRDtrapConfig::Instance();
 192     fSimParam      = AliTRDSimParam::Instance();
 193     fCommonParam   = AliTRDCommonParam::Instance();
 194     fCal           = AliTRDcalibDB::Instance();
 195     fGeo           = new AliTRDgeometry();
 196   }
 197
 198   fDetector      = det;
 199   fRobPos        = robPos;
 200   fMcmPos        = mcmPos;
 201   fNADC          = fFeeParam->GetNadcMcm();
 202   fNTimeBin      = fCal->GetNumberOfTimeBins();
 203   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
 204   fMaxTracklets  = fFeeParam->GetMaxNrOfTracklets();
 205
 206   if (!fInitialized) {
 207     fADCR    = new Int_t *[fNADC];
 208     fADCF    = new Int_t *[fNADC];
 209     fZSM     = new Int_t *[fNADC];
 210     fZSM1Dim = new Int_t  [fNADC];
 211     fGainCounterA = new UInt_t[fNADC];
 212     fGainCounterB = new UInt_t[fNADC];
 213     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 214       fADCR[iadc] = new Int_t[fNTimeBin];
 215       fADCF[iadc] = new Int_t[fNTimeBin];
 216       fZSM [iadc] = new Int_t[fNTimeBin];
 217     }
 218
 219     // filter registers
 220     fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter
 221     fTailAmplLong = new UShort_t[fNADC];
 222     fTailAmplShort = new UShort_t[fNADC];
 223
 224     // tracklet calculation
 225     fFitReg = new FitReg_t[fNADC];
 226     fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets);
 227
 228     fMCMT = new UInt_t[fMaxTracklets];
 229   }
 230
 231   fInitialized = kTRUE;
 232
 233   Reset();
 234 }
 235
 236 void AliTRDmcmSim::Reset()
 237 {
 238   // Resets the data values and internal filter registers
 239   // by re-initialising them
 240
 241   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 242     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 243       fADCR[iadc][it] = 0;
 244       fADCF[iadc][it] = 0;
 245       fZSM [iadc][it] = 1;   // Default unread = 1
 246     }
 247     fZSM1Dim[iadc] = 1;      // Default unread = 1
 248     fGainCounterA[iadc] = 0;
 249     fGainCounterB[iadc] = 0;
 250   }
 251
 252   for(Int_t i = 0; i < fMaxTracklets; i++) {
 253     fMCMT[i] = 0;
 254   }
 255
 256   FilterPedestalInit();
 257   FilterGainInit();
 258   FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
 259 }
 260
 261 void AliTRDmcmSim::SetNTimebins(Int_t ntimebins)
 262 {
 263   fNTimeBin = ntimebins;
 264   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 265     delete fADCR[iadc];
 266     delete fADCF[iadc];
 267     delete fZSM[iadc];
 268     fADCR[iadc] = new Int_t[fNTimeBin];
 269     fADCF[iadc] = new Int_t[fNTimeBin];
 270     fZSM [iadc] = new Int_t[fNTimeBin];
 271   }
 272 }
 273
 274 Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm)
 275 {
 276   // loads the ADC data as obtained from the digitsManager for the specified MCM
 277
 278   Init(det, rob, mcm);
 279
 280   if (!runloader) {
 281     AliError("No Runloader given");
 282     return kFALSE;
 283   }
 284
 285   AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
 286   if (!trdLoader) {
 287     AliError("Could not get TRDLoader");
 288     return kFALSE;
 289   }
 290
 291   Bool_t retval = kTRUE;
 292   trdLoader->LoadDigits();
 293   fDigitsManager = 0x0;
 294   AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
 295   digMgr->SetSDigits(0);
 296   digMgr->CreateArrays();
 297   digMgr->ReadDigits(trdLoader->TreeD());
 298   AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det);
 299   if (digits->HasData()) {
 300     digits->Expand();
 301
 302     if (fNTimeBin != digits->GetNtime())
 303       SetNTimebins(digits->GetNtime());
 304
 305     Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm);
 306     Int_t padcol = 0;
 307     for (Int_t ch = 0; ch < fNADC; ch++) {
 308       padcol = GetCol(ch);
 309       fZSM1Dim[ch] = 1;
 310       if (padcol < 0) {
 311         fZSM1Dim[ch] = 0;
 312         for (Int_t tb = 0; tb < fNTimeBin; tb++) {
 313           fADCR[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 314           fADCF[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 315         }
 316       }
 317       else {
 318         for (Int_t tb = 0; tb < fNTimeBin; tb++) {
 319           if (digits->GetData(padrow,padcol, tb) < 0) {
 320             fZSM1Dim[ch] = 0;
 321             fADCR[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 322             fADCF[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 323           }
 324           else {
 325             fADCR[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits + (fgAddBaseline << fgkAddDigits);
 326             fADCF[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits + (fgAddBaseline << fgkAddDigits);
 327           }
 328         }
 329       }
 330     }
 331   }
 332   else
 333     retval = kFALSE;
 334
 335   delete digMgr;
 336
 337   return retval;
 338 }
 339
 340 void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
 341 {
 342   // This function can be used to test the filters.
 343   // It feeds nsamples of ADC values with a gaussian distribution specified by mean and sigma.
 344   // The filter chain implemented here consists of:
 345   // Pedestal -> Gain -> Tail
 346   // With inputGain and inputTail the input to the gain and tail filter, respectively,
 347   // can be chosen where
 348   // 0: noise input
 349   // 1: pedestal output
 350   // 2: gain output
 351   // The input has to be chosen from a stage before.
 352   // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the
 353   // same way as in normal simulation.
 354   // The functions produces four histograms with the values at the different stages.
 355
 356   TH1F *h   = new TH1F("noise", "Gaussian Noise;sample;ADC count",
 357                        nsamples, 0, nsamples);
 358   TH1F *hfp = new TH1F("pedf", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples);
 359   TH1F *hfg = new TH1F("pedg", "Pedestal #rightarrow Gain;sample;ADC count", nsamples, 0, nsamples);
 360   TH1F *hft = new TH1F("pedt", "Gain #rightarrow Tail;sample;ADC count", nsamples, 0, nsamples);
 361   h->SetStats(kFALSE);
 362   hfp->SetStats(kFALSE);
 363   hfg->SetStats(kFALSE);
 364   hft->SetStats(kFALSE);
 365
 366   Int_t value;  // ADC count with noise (10 bit)
 367   Int_t valuep; // pedestal filter output (12 bit)
 368   Int_t valueg; // gain filter output (12 bit)
 369   Int_t valuet; // tail filter value (12 bit)
 370
 371   for (Int_t i = 0; i < nsamples; i++) {
 372     value = (Int_t) gRandom->Gaus(mean, sigma);  // generate noise with gaussian distribution
 373     h->SetBinContent(i, value);
 374
 375     valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
 376
 377     if (inputGain == 0)
 378       valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
 379     else
 380       valueg = FilterGainNextSample(1, valuep);
 381
 382     if (inputTail == 0)
 383       valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
 384     else if (inputTail == 1)
 385       valuet = FilterTailNextSample(1, valuep);
 386     else
 387       valuet = FilterTailNextSample(1, valueg);
 388
 389     hfp->SetBinContent(i, valuep >> 2);
 390     hfg->SetBinContent(i, valueg >> 2);
 391     hft->SetBinContent(i, valuet >> 2);
 392   }
 393
 394   TCanvas *c = new TCanvas;
 395   c->Divide(2,2);
 396   c->cd(1);
 397   h->Draw();
 398   c->cd(2);
 399   hfp->Draw();
 400   c->cd(3);
 401   hfg->Draw();
 402   c->cd(4);
 403   hft->Draw();
 404 }
 405
 406 Bool_t AliTRDmcmSim::CheckInitialized()
 407 {
 408   //
 409   // Check whether object is initialized
 410   //
 411
 412   if( ! fInitialized ) {
 413     AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
 414   }
 415   return fInitialized;
 416 }
 417
 418 void AliTRDmcmSim::Print(Option_t* const option) const
 419 {
 420   // Prints the data stored and/or calculated for this MCM.
 421   // The output is controlled by option which can be a sequence of any of
 422   // the following characters:
 423   // R - prints raw ADC data
 424   // F - prints filtered data
 425   // H - prints detected hits
 426   // T - prints found tracklets
 427   // The later stages are only useful when the corresponding calculations
 428   // have been performed.
 429
 430   printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
 431
 432   TString opt = option;
 433   if (opt.Contains("R")) {
 434     printf("Raw ADC data (10 bit):\n");
 435     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 436       for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
 437         printf("%5i", fADCR[iChannel][iTimeBin] >> fgkAddDigits);
 438       }
 439       printf("\n");
 440     }
 441   }
 442
 443   if (opt.Contains("F")) {
 444     printf("Filtered data (12 bit):\n");
 445     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 446       for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
 447         printf("%5i", fADCF[iChannel][iTimeBin]);
 448       }
 449       printf("\n");
 450     }
 451   }
 452
 453   if (opt.Contains("H")) {
 454     printf("Found %i hits:\n", fNHits);
 455     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
 456       printf("Hit %3i in timebin %2i, ADC %2i has charge %3i and position %3i\n",
 457              iHit,  fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
 458     }
 459   }
 460
 461   if (opt.Contains("T")) {
 462     printf("Tracklets:\n");
 463     for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntriesFast(); iTrkl++) {
 464       printf("tracklet %i: 0x%08x\n", iTrkl, ((AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl])->GetTrackletWord());
 465     }
 466   }
 467 }
 468
 469 void AliTRDmcmSim::Draw(Option_t* const option)
 470 {
 471   // Plots the data stored in a 2-dim. timebin vs. ADC channel plot.
 472   // The option selects what data is plotted and can be a sequence of
 473   // the following characters:
 474   // R - plot raw data (default)
 475   // F - plot filtered data (meaningless if R is specified)
 476   // In addition to the ADC values:
 477   // H - plot hits
 478   // T - plot tracklets
 479
 480   TString opt = option;
 481
 482   TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \
 483                                         fMcmPos, fRobPos, fDetector), \
 484                         fNADC, -0.5, fNADC-.5, fNTimeBin, -.5, fNTimeBin-.5);
 485   hist->GetXaxis()->SetTitle("ADC Channel");
 486   hist->GetYaxis()->SetTitle("Timebin");
 487   hist->SetStats(kFALSE);
 488
 489   if (opt.Contains("R")) {
 490     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 491       for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 492         hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits);
 493       }
 494     }
 495   }
 496   else {
 497     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 498       for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 499         hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
 500       }
 501     }
 502   }
 503   hist->Draw("colz");
 504
 505   if (opt.Contains("H")) {
 506     TGraph *grHits = new TGraph();
 507     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
 508       grHits->SetPoint(iHit,
 509                        fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256.,
 510                        fHits[iHit].fTimebin);
 511     }
 512     grHits->Draw("*");
 513   }
 514
 515   if (opt.Contains("T")) {
 516     TLine *trklLines = new TLine[4];
 517     for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
 518       AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
 519       AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
 520       Float_t offset = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19)) + 19 * pp->GetWidthIPad();
 521       trklLines[iTrkl].SetX1((offset -  trkl->GetY()) / pp->GetWidthIPad());
 522       trklLines[iTrkl].SetY1(0);
 523       trklLines[iTrkl].SetX2((offset - (trkl->GetY() + ((Float_t) trkl->GetdY())*140e-4)) / pp->GetWidthIPad());
 524       trklLines[iTrkl].SetY2(fNTimeBin - 1);
 525       trklLines[iTrkl].SetLineColor(2);
 526       trklLines[iTrkl].SetLineWidth(2);
 527       printf("Tracklet %i: y = %f, dy = %f, offset = %f\n", iTrkl, trkl->GetY(), (trkl->GetdY() * 140e-4), offset);
 528       trklLines[iTrkl].Draw();
 529     }
 530   }
 531 }
 532
 533 void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
 534 {
 535   //
 536   // Store ADC data into array of raw data
 537   //
 538
 539   if( !CheckInitialized() ) return;
 540
 541   if( iadc < 0 || iadc >= fNADC ) {
 542                         //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 543     return;
 544   }
 545
 546   for( Int_t it = 0 ;  it < fNTimeBin ; it++ ) {
 547     fADCR[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
 548     fADCF[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
 549   }
 550 }
 551
 552 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
 553 {
 554   //
 555   // Store ADC data into array of raw data
 556   //
 557
 558   if( !CheckInitialized() ) return;
 559
 560   if( iadc < 0 || iadc >= fNADC ) {
 561     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 562     return;
 563   }
 564
 565   fADCR[iadc][it] = adc << fgkAddDigits;
 566   fADCF[iadc][it] = adc << fgkAddDigits;
 567 }
 568
 569 void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
 570 {
 571   // Set the ADC data from an AliTRDarrayADC
 572
 573   if (!fInitialized) {
 574     AliError("Called uninitialized! Nothing done!");
 575     return;
 576   }
 577
 578   fDigitsManager = digitsManager;
 579
 580   if (fNTimeBin != adcArray->GetNtime())
 581     SetNTimebins(adcArray->GetNtime());
 582
 583   Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
 584
 585 //  Int_t firstAdc = 0;
 586 //  Int_t lastAdc = fNADC-1;
 587 //
 588 //  while (GetCol(firstAdc) < 0) {
 589 //    for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 590 //      fADCR[firstAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 591 //      fADCF[firstAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 592 //    }
 593 //    firstAdc++;
 594 //  }
 595 //
 596 //  while (GetCol(lastAdc) < 0) {
 597 //    for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 598 //      fADCR[lastAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 599 //      fADCF[lastAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 600 //    }
 601 //    lastAdc--;
 602 //  }
 603
 604   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 605     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 606       Int_t value = adcArray->GetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin);
 607       if (value < 0 || (20-iAdc + offset < 1) || (20-iAdc + offset > 165)) {
 608         fADCR[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 609         fADCF[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 610       }
 611       else {
 612         fADCR[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits + (fgAddBaseline << fgkAddDigits);
 613         fADCF[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits + (fgAddBaseline << fgkAddDigits);
 614       }
 615     }
 616   }
 617 }
 618
 619 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
 620 {
 621   //
 622   // Store ADC data into array of raw data
 623   //
 624
 625   if( !CheckInitialized() ) return;
 626
 627   if( iadc < 0 || iadc >= fNADC ) {
 628     return;
 629   }
 630
 631   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 632     fADCR[iadc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 633     fADCF[iadc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
 634   }
 635 }
 636
 637 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
 638 {
 639   //
 640   // Return column id of the pad for the given ADC channel
 641   //
 642
 643   if( !CheckInitialized() )
 644     return -1;
 645
 646   Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
 647   if (col < 0 || col >= fFeeParam->GetNcol())
 648     return -1;
 649   else
 650     return col;
 651 }
 652
 653 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
 654 {
 655   //
 656   // Produce raw data stream from this MCM and put in buf
 657   // Returns number of words filled, or negative value
 658   // with -1 * number of overflowed words
 659   //
 660
 661   UInt_t  x;
 662   Int_t   nw  = 0;  // Number of written words
 663   Int_t   of  = 0;  // Number of overflowed words
 664   Int_t   rawVer   = fFeeParam->GetRAWversion();
 665   Int_t **adc;
 666   Int_t   nActiveADC = 0;       // number of activated ADC bits in a word
 667
 668   if( !CheckInitialized() ) return 0;
 669
 670   if( fFeeParam->GetRAWstoreRaw() ) {
 671     adc = fADCR;
 672   } else {
 673     adc = fADCF;
 674   }
 675
 676   // Produce MCM header
 677   x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
 678
 679   if (nw < maxSize) {
 680     buf[nw++] = x;
 681     //printf("\nMCM header: %X ",x);
 682   }
 683   else {
 684     of++;
 685   }
 686
 687   // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
 688   //                            n : unused , c : ADC count, m : selected ADCs
 689   if( rawVer >= 3 ) {
 690     x = 0;
 691     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
 692       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
 693                 x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
 694                 nActiveADC++;           // number of 1 in mmm....m
 695       }
 696     }
 697         x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
 698         //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
 699
 700     if (nw < maxSize) {
 701       buf[nw++] = x;
 702       //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
 703     }
 704     else {
 705       of++;
 706     }
 707   }
 708
 709   // Produce ADC data. 3 timebins are packed into one 32 bits word
 710   // In this version, different ADC channel will NOT share the same word
 711
 712   UInt_t aa=0, a1=0, a2=0, a3=0;
 713
 714   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
 715     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
 716     aa = !(iAdc & 1) + 2;
 717     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
 718       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] >> fgkAddDigits : 0;
 719       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
 720       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
 721       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
 722       if (nw < maxSize) {
 723         buf[nw++] = x;
 724         //printf("%08X ",x);
 725       }
 726       else {
 727         of++;
 728       }
 729     }
 730   }
 731
 732   if( of != 0 ) return -of; else return nw;
 733 }
 734
 735 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
 736 {
 737   //
 738   // Produce tracklet data stream from this MCM and put in buf
 739   // Returns number of words filled, or negative value
 740   // with -1 * number of overflowed words
 741   //
 742
 743   Int_t   nw  = 0;  // Number of written words
 744   Int_t   of  = 0;  // Number of overflowed words
 745
 746   if( !CheckInitialized() ) return 0;
 747
 748   // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
 749   // fMCMT is filled continuously until no more tracklet words available
 750
 751   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
 752     if (nw < maxSize)
 753       buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
 754     else
 755       of++;
 756   }
 757
 758   if( of != 0 ) return -of; else return nw;
 759 }
 760
 761 void AliTRDmcmSim::Filter()
 762 {
 763   //
 764   // Filter the raw ADC values. The active filter stages and their
 765   // parameters are taken from AliTRDtrapConfig.
 766   // The raw data is stored separate from the filtered data. Thus,
 767   // it is possible to run the filters on a set of raw values
 768   // sequentially for parameter tuning.
 769   //
 770
 771   if( !CheckInitialized() ) {
 772     AliError("got called before initialization! Nothing done!");
 773     return;
 774   }
 775
 776   // Apply filters sequentially. Bypass is handled by filters
 777   // since counters and internal registers may be updated even
 778   // if the filter is bypassed.
 779   // The first filter takes the data from fADCR and
 780   // outputs to fADCF.
 781
 782   // Non-linearity filter not implemented.
 783   FilterPedestal();
 784   FilterGain();
 785   FilterTail();
 786   // Crosstalk filter not implemented.
 787 }
 788
 789 void AliTRDmcmSim::FilterPedestalInit()
 790 {
 791   // Initializes the pedestal filter assuming that the input has
 792   // been constant for a long time (compared to the time constant).
 793
 794 //  UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
 795   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
 796   UShort_t    shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from?
 797
 798   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
 799     fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
 800 }
 801
 802 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
 803 {
 804   // Returns the output of the pedestal filter given the input value.
 805   // The output depends on the internal registers and, thus, the
 806   // history of the filter.
 807
 808   UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
 809   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
 810   UShort_t    fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low
 811   UShort_t    shifts[4] = {11, 14, 17, 21}; //??? where to come from
 812
 813   UShort_t accumulatorShifted;
 814   Int_t correction;
 815   UShort_t inpAdd;
 816
 817   inpAdd = value + fpnp;
 818
 819   if (fpby == 0) //??? before or after update of accumulator
 820     return value;
 821
 822   accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF;   // 10 bits
 823   if (timebin == 0) // the accumulator is disabled in the drift time
 824   {
 825     correction = (value & 0x3FF) - accumulatorShifted;
 826     fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF;             // 31 bits
 827   }
 828
 829   if (inpAdd <= accumulatorShifted)
 830     return 0;
 831   else
 832   {
 833     inpAdd = inpAdd - accumulatorShifted;
 834     if (inpAdd > 0xFFF)
 835       return 0xFFF;
 836     else
 837       return inpAdd;
 838   }
 839 }
 840
 841 void AliTRDmcmSim::FilterPedestal()
 842 {
 843   //
 844   // Apply pedestal filter
 845   //
 846   // As the first filter in the chain it reads data from fADCR
 847   // and outputs to fADCF.
 848   // It has only an effect if previous samples have been fed to
 849   // find the pedestal. Currently, the simulation assumes that
 850   // the input has been stable for a sufficiently long time.
 851
 852   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 853     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 854       fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
 855     }
 856   }
 857 }
 858
 859 void AliTRDmcmSim::FilterGainInit()
 860 {
 861   // Initializes the gain filter. In this case, only threshold
 862   // counters are reset.
 863
 864   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 865     // these are counters which in hardware continue
 866     // until maximum or reset
 867     fGainCounterA[iAdc] = 0;
 868     fGainCounterB[iAdc] = 0;
 869   }
 870 }
 871
 872 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
 873 {
 874   // Apply the gain filter to the given value.
 875   // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
 876   // The output depends on the internal registers and, thus, the
 877   // history of the filter.
 878
 879   UShort_t    fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
 880   UShort_t    fgf  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
 881   UShort_t    fga  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
 882   UShort_t    fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
 883   UShort_t    fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
 884
 885   UInt_t tmp;
 886
 887   value &= 0xFFF;
 888   tmp = (value * fgf) >> 11;
 889   if (tmp > 0xFFF) tmp = 0xFFF;
 890
 891   if (fgby == 1)
 892     value = AddUintClipping(tmp, fga, 12);
 893
 894   // Update threshold counters
 895   // not really useful as they are cleared with every new event
 896   if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
 897   {
 898     if (value >= fgtb)
 899       fGainCounterB[adc]++;
 900     else if (value >= fgta)
 901       fGainCounterA[adc]++;
 902   }
 903
 904   return value;
 905 }
 906
 907 void AliTRDmcmSim::FilterGain()
 908 {
 909   // Read data from fADCF and apply gain filter.
 910
 911   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 912     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
 913         fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
 914     }
 915   }
 916 }
 917
 918 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
 919 {
 920   // Initializes the tail filter assuming that the input has
 921   // been at the baseline value (configured by FTFP) for a
 922   // sufficiently long time.
 923
 924   // exponents and weight calculated from configuration
 925   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
 926   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
 927   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
 928
 929   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
 930   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
 931   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
 932   Float_t qup, qdn;
 933   qup = (1 - lambdaL) * (1 - lambdaS);
 934   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
 935   Float_t kdc = qup/qdn;
 936
 937   Float_t kt, ql, qs;
 938   UShort_t aout;
 939
 940   kt = kdc * baseline;
 941   aout = baseline - (UShort_t) kt;
 942   ql = lambdaL * (1 - lambdaS) *      alphaL;
 943   qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
 944
 945   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
 946     fTailAmplLong[iAdc]  = (UShort_t) (aout * ql / (ql + qs));
 947     fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
 948   }
 949 }
 950
 951 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
 952 {
 953   // Returns the output of the tail filter for the given input value.
 954   // The output depends on the internal registers and, thus, the
 955   // history of the filter.
 956
 957   // exponents and weight calculated from configuration
 958   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
 959   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
 960   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
 961
 962   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
 963   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
 964   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
 965   Float_t qup, qdn;
 966   qup = (1 - lambdaL) * (1 - lambdaS);
 967   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
 968 //  Float_t kdc = qup/qdn;
 969
 970   UInt_t aDiff;
 971   UInt_t alInpv;
 972   UShort_t aQ;
 973   UInt_t tmp;
 974
 975   UShort_t inpVolt = value & 0xFFF;    // 12 bits
 976
 977   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
 978     return value;
 979   else
 980   {
 981     // add the present generator outputs
 982     aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
 983
 984     // calculate the difference between the input the generated signal
 985     if (inpVolt > aQ)
 986       aDiff = inpVolt - aQ;
 987     else
 988       aDiff = 0;
 989
 990     // the inputs to the two generators, weighted
 991     alInpv = (aDiff * alphaLong) >> 11;
 992
 993     // the new values of the registers, used next time
 994     // long component
 995     tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
 996     tmp =  (tmp * lambdaLong) >> 11;
 997     fTailAmplLong[adc] = tmp & 0xFFF;
 998     // short component
 999     tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
1000     tmp =  (tmp * lambdaShort) >> 11;
1001     fTailAmplShort[adc] = tmp & 0xFFF;
1002
1003     // the output of the filter
1004     return aDiff;
1005   }
1006 }
1007
1008 void AliTRDmcmSim::FilterTail()
1009 {
1010   // Apply tail filter
1011
1012   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1013     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
1014       fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
1015     }
1016   }
1017 }
1018
1019 void AliTRDmcmSim::ZSMapping()
1020 {
1021   //
1022   // Zero Suppression Mapping implemented in TRAP chip
1023   //
1024   // See detail TRAP manual "Data Indication" section:
1025   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
1026   //
1027
1028   //??? values should come from TRAPconfig
1029   Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4  (Tis=4)
1030   Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40)
1031   Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0
1032                                                                  // (lookup table accept (I2,I1,I0)=(111)
1033                                                                  // or (110) or (101) or (100))
1034   Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
1035   Int_t ep   = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here
1036
1037   Int_t **adc = fADCF;
1038
1039   if( !CheckInitialized() ) {
1040     AliError("got called uninitialized! Nothing done!");
1041     return;
1042   }
1043
1044   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1045     for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
1046
1047       // Get ADC data currently in filter buffer
1048       Int_t ap = adc[iadc-1][it] - ep; // previous
1049       Int_t ac = adc[iadc  ][it] - ep; // current
1050       Int_t an = adc[iadc+1][it] - ep; // next
1051
1052       // evaluate three conditions
1053       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
1054       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
1055       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
1056
1057       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
1058       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
1059                                          // and d=1 means false according to TRAP manual)
1060
1061       fZSM[iadc][it] &= d;
1062       if( eBIN == 0 ) {  // turn on neighboring ADCs
1063         fZSM[iadc-1][it] &= d;
1064         fZSM[iadc+1][it] &= d;
1065       }
1066     }
1067   }
1068
1069   // do 1 dim projection
1070   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1071     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1072       fZSM1Dim[iadc] &= fZSM[iadc][it];
1073     }
1074   }
1075 }
1076
1077 void AliTRDmcmSim::DumpData( const char * const f, const char * const target )
1078 {
1079   //
1080   // Dump data stored (for debugging).
1081   // target should contain one or multiple of the following characters
1082   //   R   for raw data
1083   //   F   for filtered data
1084   //   Z   for zero suppression map
1085   //   S   Raw dat astream
1086   // other characters are simply ignored
1087   //
1088
1089   UInt_t tempbuf[1024];
1090
1091   if( !CheckInitialized() ) return;
1092
1093   std::ofstream of( f, std::ios::out | std::ios::app );
1094   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
1095              fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector),
1096              fGeo->GetSector(fDetector), fRobPos, fMcmPos );
1097
1098   for( Int_t t=0 ; target[t] != 0 ; t++ ) {
1099     switch( target[t] ) {
1100     case 'R' :
1101     case 'r' :
1102       of << Form("fADCR (raw ADC data)\n");
1103       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1104         of << Form("  ADC %02d: ", iadc);
1105         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1106           of << Form("% 4d",  fADCR[iadc][it]);
1107         }
1108         of << Form("\n");
1109       }
1110       break;
1111     case 'F' :
1112     case 'f' :
1113       of << Form("fADCF (filtered ADC data)\n");
1114       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1115         of << Form("  ADC %02d: ", iadc);
1116         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1117           of << Form("% 4d",  fADCF[iadc][it]);
1118         }
1119         of << Form("\n");
1120       }
1121       break;
1122     case 'Z' :
1123     case 'z' :
1124       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
1125       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1126         of << Form("  ADC %02d: ", iadc);
1127         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
1128         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1129           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
1130         }
1131         of << Form("\n");
1132       }
1133       break;
1134     case 'S' :
1135     case 's' :
1136       Int_t s = ProduceRawStream( tempbuf, 1024 );
1137       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
1138       of << Form("  address  data\n");
1139       for( Int_t i = 0 ; i < s ; i++ ) {
1140         of << Form("  %04x     %08x\n", i, tempbuf[i]);
1141       }
1142     }
1143   }
1144 }
1145
1146 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
1147 {
1148   // Add the given hit to the fit register which is lateron used for
1149   // the tracklet calculation.
1150   // In addition to the fit sums in the fit register MC information
1151   // is stored.
1152
1153   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) &&
1154       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1155     fFitReg[adc].fQ0 += qtot;
1156
1157   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1158       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1159     fFitReg[adc].fQ1 += qtot;
1160
1161   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1162       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1163   {
1164     fFitReg[adc].fSumX  += timebin;
1165     fFitReg[adc].fSumX2 += timebin*timebin;
1166     fFitReg[adc].fNhits++;
1167     fFitReg[adc].fSumY  += ypos;
1168     fFitReg[adc].fSumY2 += ypos*ypos;
1169     fFitReg[adc].fSumXY += timebin*ypos;
1170   }
1171
1172   // register hits (MC info)
1173   fHits[fNHits].fChannel = adc;
1174   fHits[fNHits].fQtot = qtot;
1175   fHits[fNHits].fYpos = ypos;
1176   fHits[fNHits].fTimebin = timebin;
1177   fHits[fNHits].fLabel = label;
1178   fNHits++;
1179 }
1180
1181 void AliTRDmcmSim::CalcFitreg()
1182 {
1183   // Preprocessing.
1184   // Detect the hits and fill the fit registers.
1185   // Requires 12-bit data from fADCF which means Filter()
1186   // has to be called before even if all filters are bypassed.
1187
1188   //???
1189   // TRAP parameters:
1190   const UShort_t lutPos[128] = {   // move later to some other file
1191     0,  1,  1,  2,  2,  3,  3,  4,  4,  5,  5,  6,  6,  7,  7,  8,  8,  9,  9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
1192     16, 16, 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26,
1193     27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, 27, 27, 26,
1194     26, 26, 26, 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 18, 18, 17, 17, 16, 15, 14, 13, 12, 11, 10,  9,  8,  7,  7};
1195
1196   //??? to be clarified:
1197   UInt_t adcMask = 0xffffffff;
1198
1199   UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
1200   Short_t ypos, fromLeft, fromRight, found;
1201   UShort_t qTotal[19]; // the last is dummy
1202   UShort_t marked[6], qMarked[6], worse1, worse2;
1203
1204   timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
1205   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)
1206       < timebin1)
1207     timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1208   timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
1209   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)
1210       > timebin2)
1211     timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1212
1213   // reset the fit registers
1214   fNHits = 0;
1215   for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
1216   {
1217     fFitReg[adcch].fNhits = 0;
1218     fFitReg[adcch].fQ0    = 0;
1219     fFitReg[adcch].fQ1    = 0;
1220     fFitReg[adcch].fSumX  = 0;
1221     fFitReg[adcch].fSumY  = 0;
1222     fFitReg[adcch].fSumX2 = 0;
1223     fFitReg[adcch].fSumY2 = 0;
1224     fFitReg[adcch].fSumXY = 0;
1225   }
1226
1227   for (timebin = timebin1; timebin < timebin2; timebin++)
1228   {
1229     // first find the hit candidates and store the total cluster charge in qTotal array
1230     // in case of not hit store 0 there.
1231     for (adcch = 0; adcch < fNADC-2; adcch++) {
1232       if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
1233       {
1234         adcLeft  = fADCF[adcch  ][timebin];
1235         adcCentral  = fADCF[adcch+1][timebin];
1236         adcRight = fADCF[adcch+2][timebin];
1237         if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1)
1238           hitQual = ( (adcLeft * adcRight) <
1239                        (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
1240         else
1241           hitQual = 1;
1242         // The accumulated charge is with the pedestal!!!
1243         qtotTemp = adcLeft + adcCentral + adcRight;
1244         if ( (hitQual) &&
1245              (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1246              (adcLeft <= adcCentral) &&
1247              (adcCentral > adcRight) )
1248           qTotal[adcch] = qtotTemp;
1249         else
1250           qTotal[adcch] = 0;
1251       }
1252       else
1253         qTotal[adcch] = 0; //jkl
1254       AliDebug(10,Form("ch %2d   qTotal %5d",adcch, qTotal[adcch]));
1255     }
1256
1257     fromLeft = -1;
1258     adcch = 0;
1259     found = 0;
1260     marked[4] = 19; // invalid channel
1261     marked[5] = 19; // invalid channel
1262     qTotal[19] = 0;
1263     while ((adcch < 16) && (found < 3))
1264     {
1265       if (qTotal[adcch] > 0)
1266       {
1267         fromLeft = adcch;
1268         marked[2*found+1]=adcch;
1269         found++;
1270       }
1271       adcch++;
1272     }
1273
1274     fromRight = -1;
1275     adcch = 18;
1276     found = 0;
1277     while ((adcch > 2) && (found < 3))
1278     {
1279       if (qTotal[adcch] > 0)
1280       {
1281         marked[2*found]=adcch;
1282         found++;
1283         fromRight = adcch;
1284       }
1285       adcch--;
1286     }
1287
1288     AliDebug(10,Form("Fromleft=%d, Fromright=%d",fromLeft, fromRight));
1289     // here mask the hit candidates in the middle, if any
1290     if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1291       for (adcch = fromLeft+1; adcch < fromRight; adcch++)
1292         qTotal[adcch] = 0;
1293
1294     found = 0;
1295     for (adcch = 0; adcch < 19; adcch++)
1296       if (qTotal[adcch] > 0) found++;
1297     // NOT READY
1298
1299     if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1300     {
1301       if (marked[4] == marked[5]) marked[5] = 19;
1302       for (found=0; found<6; found++)
1303       {
1304         qMarked[found] = qTotal[marked[found]] >> 4;
1305         AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found]));
1306       }
1307
1308       Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1309                     qMarked[0],
1310                     qMarked[3],
1311                     qMarked[4],
1312                     qMarked[1],
1313                     qMarked[2],
1314                     qMarked[5],
1315                     &worse1, &worse2);
1316       // Now mask the two channels with the smallest charge
1317       if (worse1 < 19)
1318       {
1319         qTotal[worse1] = 0;
1320         AliDebug(10,Form("Kill ch %d\n",worse1));
1321       }
1322       if (worse2 < 19)
1323       {
1324         qTotal[worse2] = 0;
1325         AliDebug(10,Form("Kill ch %d\n",worse2));
1326       }
1327     }
1328
1329     for (adcch = 0; adcch < 19; adcch++) {
1330       if (qTotal[adcch] > 0) // the channel is marked for processing
1331       {
1332         adcLeft  = fADCF[adcch  ][timebin];
1333         adcCentral  = fADCF[adcch+1][timebin];
1334         adcRight = fADCF[adcch+2][timebin];
1335         // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1336         // subtract the pedestal TPFP, clipping instead of wrapping
1337
1338         Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1339         AliDebug(10, Form("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1340                timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP,
1341                fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)));
1342
1343         if (adcLeft  < regTPFP) adcLeft  = 0; else adcLeft  -= regTPFP;
1344         if (adcCentral  < regTPFP) adcCentral  = 0; else adcCentral  -= regTPFP;
1345         if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
1346
1347         // Calculate the center of gravity
1348         // checking for adcCentral != 0 (in case of "bad" configuration)
1349         if (adcCentral == 0)
1350           continue;
1351         ypos = 128*(adcLeft - adcRight) / adcCentral;
1352         if (ypos < 0) ypos = -ypos;
1353         // make the correction using the LUT
1354         ypos = ypos + lutPos[ypos & 0x7F];
1355         if (adcLeft > adcRight) ypos = -ypos;
1356
1357         // label calculation
1358         Int_t mcLabel = -1;
1359         if (fDigitsManager) {
1360           Int_t label[9] = { 0 }; // up to 9 different labels possible
1361           Int_t count[9] = { 0 };
1362           Int_t maxIdx = -1;
1363           Int_t maxCount = 0;
1364           Int_t nLabels = 0;
1365           Int_t padcol[3];
1366           padcol[0] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch);
1367           padcol[1] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+1);
1368           padcol[2] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+2);
1369           Int_t padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1370           for (Int_t iDict = 0; iDict < 3; iDict++) {
1371             if (!fDigitsManager->UsesDictionaries() || fDigitsManager->GetDictionary(fDetector, iDict) == 0) {
1372               AliError("Cannot get dictionary");
1373               continue;
1374             }
1375             AliTRDarrayDictionary *dict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict);
1376             if (dict->GetDim() == 0) {
1377               AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
1378               continue;
1379             }
1380             dict->Expand();
1381             for (Int_t iPad = 0; iPad < 3; iPad++) {
1382               if (padcol[iPad] < 0)
1383                 continue;
1384               Int_t currLabel = dict->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector);
1385               AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin));
1386               for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1387                 if (currLabel == label[iLabel]) {
1388                   count[iLabel]++;
1389                   if (count[iLabel] > maxCount) {
1390                     maxCount = count[iLabel];
1391                     maxIdx = iLabel;
1392                   }
1393                   currLabel = 0;
1394                   break;
1395                 }
1396               }
1397               if (currLabel > 0) {
1398                 label[nLabels++] = currLabel;
1399               }
1400             }
1401           }
1402           if (maxIdx >= 0)
1403             mcLabel = label[maxIdx];
1404         }
1405
1406         // add the hit to the fitregister
1407         AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel);
1408       }
1409     }
1410   }
1411 }
1412
1413 void AliTRDmcmSim::TrackletSelection()
1414 {
1415   // Select up to 4 tracklet candidates from the fit registers
1416   // and assign them to the CPUs.
1417
1418   UShort_t adcIdx, i, j, ntracks, tmp;
1419   UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
1420
1421   ntracks = 0;
1422   for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1423     if ( (fFitReg[adcIdx].fNhits
1424           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
1425          (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
1426           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1427     {
1428       trackletCand[ntracks][0] = adcIdx;
1429       trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1430       AliDebug(10,Form("%d  %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]));
1431       ntracks++;
1432     };
1433
1434   for (i=0; i<ntracks;i++)
1435     AliDebug(10,Form("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]));
1436
1437   if (ntracks > 4)
1438   {
1439     // primitive sorting according to the number of hits
1440     for (j = 0; j < (ntracks-1); j++)
1441     {
1442       for (i = j+1; i < ntracks; i++)
1443       {
1444         if ( (trackletCand[j][1]  < trackletCand[i][1]) ||
1445              ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
1446         {
1447           // swap j & i
1448           tmp = trackletCand[j][1];
1449           trackletCand[j][1] = trackletCand[i][1];
1450           trackletCand[i][1] = tmp;
1451           tmp = trackletCand[j][0];
1452           trackletCand[j][0] = trackletCand[i][0];
1453           trackletCand[i][0] = tmp;
1454         }
1455       }
1456     }
1457     ntracks = 4; // cut the rest, 4 is the max
1458   }
1459   // else is not necessary to sort
1460
1461   // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1462   for (j = 0; j < (ntracks-1); j++)
1463   {
1464     for (i = j+1; i < ntracks; i++)
1465     {
1466       if (trackletCand[j][0] < trackletCand[i][0])
1467       {
1468         // swap j & i
1469         tmp = trackletCand[j][1];
1470         trackletCand[j][1] = trackletCand[i][1];
1471         trackletCand[i][1] = tmp;
1472         tmp = trackletCand[j][0];
1473         trackletCand[j][0] = trackletCand[i][0];
1474         trackletCand[i][0] = tmp;
1475       }
1476     }
1477   }
1478   for (i = 0; i < ntracks; i++)  // CPUs with tracklets.
1479     fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
1480   for (i = ntracks; i < 4; i++)  // CPUs without tracklets
1481     fFitPtr[i] = 31;            // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1482   AliDebug(10,Form("found %i tracklet candidates\n", ntracks));
1483   for (i = 0; i < 4; i++)
1484     AliDebug(10,Form("fitPtr[%i]: %i\n", i, fFitPtr[i]));
1485 }
1486
1487 void AliTRDmcmSim::FitTracklet()
1488 {
1489   // Perform the actual tracklet fit based on the fit sums
1490   // which have been filled in the fit registers.
1491
1492   // parameters in fitred.asm (fit program)
1493   Int_t decPlaces = 5;
1494   Int_t rndAdd = 0;
1495   if (decPlaces >  1)
1496     rndAdd = (1 << (decPlaces-1)) + 1;
1497   else if (decPlaces == 1)
1498     rndAdd = 1;
1499   Int_t ndriftDp = 5;  // decimal places for drift time
1500   Long64_t shift = ((Long64_t) 1 << 32);
1501
1502
1503   // calculated in fitred.asm
1504   Int_t padrow = ((fRobPos >> 1) << 2) | (fMcmPos >> 2);
1505   Int_t yoffs = (((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) -
1506     ((18*4*2 - 18*2 - 1) << 7);
1507   yoffs = yoffs << decPlaces; // holds position of ADC channel 1
1508   Int_t layer = fDetector % 6;
1509   UInt_t scaleY = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 160.0e-4) * shift);
1510   UInt_t scaleD = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 140.0e-4) * shift);
1511   // previously taken from geometry:
1512   // UInt_t scaleYold = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
1513   // UInt_t scaleDold = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
1514
1515
1516   // should come from trapConfig (DMEM)
1517   AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
1518   Float_t scaleSlope = (256 / pp->GetWidthIPad()) * (1 << decPlaces); // only used for calculation of corrections and cut
1519   Int_t ndrift   = 20 << ndriftDp; //??? value in simulation?
1520   Int_t deflCorr = (Int_t) (TMath::Tan(fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector))) * fGeo->CdrHght() * scaleSlope); // -370;
1521   Int_t tiltCorr = (Int_t) (pp->GetRowPos(padrow) / fGeo->GetTime0(fDetector % 6) * fGeo->CdrHght() * scaleSlope *
1522                             TMath::Tan(pp->GetTiltingAngle() / 180. * TMath::Pi()));
1523 //  printf("vdrift av.: %f\n", fCal->GetVdriftAverage(fDetector));
1524 //  printf("chamber height: %f\n", fGeo->CdrHght());
1525 //  printf("omega tau: %f\n", fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector)));
1526 //  printf("deflection correction: %i\n", deflCorr);
1527   Float_t ptcut = 2.3;
1528   AliMagF* fld = (AliMagF *) TGeoGlobalMagField::Instance()->GetField();
1529   Double_t bz = 0;
1530   if (fld) {
1531     bz       = 0.1 * fld->SolenoidField();   // kGauss -> Tesla
1532   }
1533 //  printf("Bz: %f\n", bz);
1534   Float_t x0 = fGeo->GetTime0(fDetector % 6);
1535   Float_t y0 = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10));
1536   Float_t alphaMax = TMath::ASin( (TMath::Sqrt(TMath::Power(x0/100., 2) + TMath::Power(y0/100., 2)) *
1537                                    0.3 * TMath::Abs(bz) ) / (2 * ptcut));
1538 //  printf("alpha max: %f\n", alphaMax * 180/TMath::Pi());
1539   Int_t minslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) + alphaMax) / 140.e-4);
1540   Int_t maxslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) - alphaMax) / 140.e-4);
1541
1542
1543   // local variables for calculation
1544   Long64_t mult, temp, denom; //???
1545   UInt_t q0, q1, qTotal;          // charges in the two windows and total charge
1546   UShort_t nHits;                 // number of hits
1547   Int_t slope, offset;            // slope and offset of the tracklet
1548   Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1549   //int32_t SumY2;                // not used in the current TRAP program
1550   FitReg_t *fit0, *fit1;          // pointers to relevant fit registers
1551
1552 //  const uint32_t OneDivN[32] = {  // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1553 //      0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1554 //      0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1555 //      0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1556 //      0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1557
1558   for (Int_t cpu = 0; cpu < 4; cpu++) {
1559     if (fFitPtr[cpu] == 31)
1560     {
1561       fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1562     }
1563     else
1564     {
1565       fit0 = &fFitReg[fFitPtr[cpu]  ];
1566       fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1567
1568       mult = 1;
1569       mult = mult << (32 + decPlaces);
1570       mult = -mult;
1571
1572       // Merging
1573       nHits   = fit0->fNhits + fit1->fNhits; // number of hits
1574       sumX    = fit0->fSumX  + fit1->fSumX;
1575       sumX2   = fit0->fSumX2 + fit1->fSumX2;
1576       denom   = nHits*sumX2 - sumX*sumX;
1577
1578       mult    = mult / denom; // exactly like in the TRAP program
1579       q0      = fit0->fQ0    + fit1->fQ0;
1580       q1      = fit0->fQ1    + fit1->fQ1;
1581       sumY    = fit0->fSumY  + fit1->fSumY  + 256*fit1->fNhits;
1582       sumXY   = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
1583
1584       slope   = nHits*sumXY - sumX * sumY;
1585       AliDebug(5, Form("slope from fitreg: %i", slope));
1586       offset  = sumX2*sumY  - sumX * sumXY;
1587       temp    = mult * slope;
1588       slope   = temp >> 32; // take the upper 32 bits
1589       slope   = -slope;
1590       temp    = mult * offset;
1591       offset  = temp >> 32; // take the upper 32 bits
1592
1593       offset = offset + yoffs;
1594       AliDebug(5, Form("slope: %i, slope * ndrift: %i, deflCorr: %i, tiltCorr: %i", slope, slope * ndrift, deflCorr, tiltCorr));
1595       slope  = ((slope * ndrift) >> ndriftDp) + deflCorr + tiltCorr;
1596       offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1597
1598       AliDebug(5, Form("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i", fDetector, fRobPos, fMcmPos, slope, minslope, maxslope));
1599       temp    = slope;
1600       temp    = temp * scaleD;
1601       slope   = (temp >> 32);
1602       AliDebug(5, Form("slope after scaling: %i", slope));
1603
1604       temp    = offset;
1605       temp    = temp * scaleY;
1606       offset  = (temp >> 32);
1607
1608       // rounding, like in the TRAP
1609       slope   = (slope  + rndAdd) >> decPlaces;
1610       AliDebug(5, Form("slope after shifting: %i", slope));
1611       offset  = (offset + rndAdd) >> decPlaces;
1612
1613       Bool_t rejected = kFALSE;
1614       if ((slope < minslope) || (slope > maxslope))
1615         rejected = kTRUE;
1616
1617       if (rejected && GetApplyCut())
1618       {
1619         fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1620       }
1621       else
1622       {
1623         if (slope > 63 || slope < -64) { // wrapping in TRAP!
1624           AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
1625           fMCMT[cpu] = 0x10001000;
1626           continue;
1627         }
1628
1629         slope   = slope  &   0x7F; // 7 bit
1630
1631         if (offset > 0xfff || offset < -0xfff)
1632           AliWarning("Overflow in offset");
1633         offset  = offset & 0x1FFF; // 13 bit
1634
1635         Float_t length = TMath::Sqrt(1 + (pp->GetRowPos(padrow) * pp->GetRowPos(padrow) +
1636                                           (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10) * pp->GetWidthIPad() *
1637                                            fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10) * pp->GetWidthIPad())) /
1638                                      (fGeo->GetTime0(fDetector % 6)*fGeo->GetTime0(fDetector % 6)));
1639
1640         //        qTotal  = (q1 / nHits) >> 1;
1641         qTotal = GetPID(q0/length/fgChargeNorm, q1/length/fgChargeNorm);
1642         if (qTotal > 0xff)
1643           AliWarning("Overflow in charge");
1644         qTotal  = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1645
1646         // assemble and store the tracklet word
1647         fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1648
1649         // calculate MC label
1650         Int_t mcLabel = -1;
1651         Int_t nHits0 = 0;
1652         Int_t nHits1 = 0;
1653         if (fDigitsManager) {
1654           Int_t label[30] = {0}; // up to 30 different labels possible
1655           Int_t count[30] = {0};
1656           Int_t maxIdx = -1;
1657           Int_t maxCount = 0;
1658           Int_t nLabels = 0;
1659           for (Int_t iHit = 0; iHit < fNHits; iHit++) {
1660             if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
1661                 (fHits[iHit].fChannel - fFitPtr[cpu] > 1))
1662               continue;
1663
1664             // counting contributing hits
1665             if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) &&
1666                 fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0))
1667               nHits0++;
1668             if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1) &&
1669                 fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
1670               nHits1++;
1671
1672             Int_t currLabel = fHits[iHit].fLabel;
1673             for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1674               if (currLabel == label[iLabel]) {
1675                 count[iLabel]++;
1676                 if (count[iLabel] > maxCount) {
1677                   maxCount = count[iLabel];
1678                   maxIdx = iLabel;
1679                 }
1680                 currLabel = 0;
1681                 break;
1682               }
1683             }
1684             if (currLabel > 0) {
1685               label[nLabels++] = currLabel;
1686             }
1687           }
1688           if (maxIdx >= 0)
1689             mcLabel = label[maxIdx];
1690         }
1691         new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1692         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
1693
1694
1695         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits(fit0->fNhits + fit1->fNhits);
1696         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits0(nHits0);
1697         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1);
1698         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
1699         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
1700       }
1701     }
1702   }
1703 }
1704
1705 Int_t AliTRDmcmSim::GetPID(Float_t q0, Float_t q1)
1706 {
1707   // get PID from accumulated charges q0 and q1
1708
1709   Int_t binQ0 = (Int_t) (q0 * fgPidNBinsQ0) + 1;
1710   Int_t binQ1 = (Int_t) (q1 * fgPidNBinsQ1) + 1;
1711   binQ0 = binQ0 >= fgPidNBinsQ0 ? fgPidNBinsQ0-1 : binQ0;
1712   binQ1 = binQ1 >= fgPidNBinsQ0 ? fgPidNBinsQ0-1 : binQ1;
1713
1714   return fgPidLut[binQ0*fgPidNBinsQ1+binQ1];
1715 }
1716
1717 void AliTRDmcmSim::SetPIDlut(Int_t *lut, Int_t nbinsq0, Int_t nbinsq1)
1718 {
1719   // set a user-defined PID LUT
1720
1721   if (fgPidLutDelete)
1722     delete [] fgPidLut;
1723
1724   fgPidLutDelete = kFALSE;
1725   fgPidLut = lut;
1726   fgPidNBinsQ0 = nbinsq0;
1727   fgPidNBinsQ1 = nbinsq1;
1728 }
1729
1730 void AliTRDmcmSim::SetPIDlut(TH2F *lut)
1731 {
1732   // set a user-defined PID LUT from a 2D histogram
1733
1734   if (fgPidLutDelete)
1735     delete [] fgPidLut;
1736
1737   fgPidNBinsQ0 = lut->GetNbinsX();
1738   fgPidNBinsQ1 = lut->GetNbinsY();
1739
1740   fgPidLut = new Int_t[fgPidNBinsQ0*fgPidNBinsQ1];
1741
1742   for (Int_t ix = 0; ix < fgPidNBinsQ0; ix++) {
1743     for (Int_t iy = 0; iy < fgPidNBinsQ1; iy++) {
1744       fgPidLut[ix*fgPidNBinsQ1 + iy] = (Int_t) (256. * lut->GetBinContent(ix, iy));
1745     }
1746   }
1747
1748   fgPidLutDelete = kTRUE;
1749 }
1750
1751 void AliTRDmcmSim::SetPIDlutDefault()
1752 {
1753   // use the default PID LUT
1754
1755   if (fgPidLutDelete )
1756     delete [] fgPidLut;
1757
1758   fgPidLutDelete = kFALSE;
1759   fgPidLut = *fgPidLutDefault;
1760   fgPidNBinsQ0 = 40;
1761   fgPidNBinsQ1 = 50;
1762 }
1763
1764 void AliTRDmcmSim::Tracklet()
1765 {
1766   // Run the tracklet calculation by calling sequentially:
1767   // CalcFitreg(); TrackletSelection(); FitTracklet()
1768   // and store the tracklets
1769
1770   if (!fInitialized) {
1771     AliError("Called uninitialized! Nothing done!");
1772     return;
1773   }
1774
1775   fTrackletArray->Delete();
1776
1777   CalcFitreg();
1778   if (fNHits == 0)
1779     return;
1780   TrackletSelection();
1781   FitTracklet();
1782 }
1783
1784 Bool_t AliTRDmcmSim::StoreTracklets()
1785 {
1786   // store the found tracklets via the loader
1787
1788   if (fTrackletArray->GetEntriesFast() == 0)
1789     return kTRUE;
1790
1791   AliRunLoader *rl = AliRunLoader::Instance();
1792   AliDataLoader *dl = 0x0;
1793   if (rl)
1794     dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1795   if (!dl) {
1796     AliError("Could not get the tracklets data loader!");
1797     return kFALSE;
1798   }
1799
1800   TTree *trackletTree = dl->Tree();
1801   if (!trackletTree) {
1802     dl->MakeTree();
1803     trackletTree = dl->Tree();
1804   }
1805
1806   AliTRDtrackletMCM *trkl = 0x0;
1807   TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1808   if (!trkbranch)
1809     trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1810
1811   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1812     trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1813     trkbranch->SetAddress(&trkl);
1814 //      printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
1815     trkbranch->Fill();
1816   }
1817   dl->WriteData("OVERWRITE");
1818
1819   return kTRUE;
1820 }
1821
1822 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1823 {
1824   // write back the processed data configured by EBSF
1825   // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1826   // zero-suppressed valued are written as -1 to digits
1827
1828   if (!fInitialized) {
1829     AliError("Called uninitialized! Nothing done!");
1830     return;
1831   }
1832
1833 //  Int_t firstAdc = 0;
1834 //  Int_t lastAdc = fNADC - 1;
1835 //
1836 //  while (GetCol(firstAdc) < 0)
1837 //    firstAdc++;
1838 //
1839 //  while (GetCol(lastAdc) < 0)
1840 //    lastAdc--;
1841
1842   Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
1843
1844   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1845   {
1846     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
1847       if (fZSM1Dim[iAdc] == 1) {
1848         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1849           digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
1850 //          printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin,
1851 //                 digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1852         }
1853       }
1854     }
1855   }
1856   else {
1857     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
1858       if (fZSM1Dim[iAdc] == 0) {
1859         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1860           digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
1861         }
1862       }
1863       else {
1864         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1865           digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
1866 //          printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1867         }
1868       }
1869     }
1870   }
1871 }
1872
1873 // help functions, to be cleaned up
1874
1875 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1876 {
1877   //
1878   // This function adds a and b (unsigned) and clips to
1879   // the specified number of bits.
1880   //
1881
1882   UInt_t sum = a + b;
1883   if (nbits < 32)
1884   {
1885     UInt_t maxv = (1 << nbits) - 1;;
1886     if (sum > maxv)
1887       sum = maxv;
1888   }
1889   else
1890   {
1891     if ((sum < a) || (sum < b))
1892       sum = 0xFFFFFFFF;
1893   }
1894   return sum;
1895 }
1896
1897 void AliTRDmcmSim::Sort2(UShort_t  idx1i, UShort_t  idx2i, \
1898                             UShort_t  val1i, UShort_t  val2i, \
1899                             UShort_t *idx1o, UShort_t *idx2o, \
1900                             UShort_t *val1o, UShort_t *val2o) const
1901 {
1902   // sorting for tracklet selection
1903
1904     if (val1i > val2i)
1905     {
1906         *idx1o = idx1i;
1907         *idx2o = idx2i;
1908         *val1o = val1i;
1909         *val2o = val2i;
1910     }
1911     else
1912     {
1913         *idx1o = idx2i;
1914         *idx2o = idx1i;
1915         *val1o = val2i;
1916         *val2o = val1i;
1917     }
1918 }
1919
1920 void AliTRDmcmSim::Sort3(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, \
1921                             UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, \
1922                             UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \
1923                             UShort_t *val1o, UShort_t *val2o, UShort_t *val3o)
1924 {
1925   // sorting for tracklet selection
1926
1927     Int_t sel;
1928
1929
1930     if (val1i > val2i) sel=4; else sel=0;
1931     if (val2i > val3i) sel=sel + 2;
1932     if (val3i > val1i) sel=sel + 1;
1933     //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1934     switch(sel)
1935     {
1936         case 6 : // 1 >  2  >  3            => 1 2 3
1937         case 0 : // 1 =  2  =  3            => 1 2 3 : in this case doesn't matter, but so is in hardware!
1938             *idx1o = idx1i;
1939             *idx2o = idx2i;
1940             *idx3o = idx3i;
1941             *val1o = val1i;
1942             *val2o = val2i;
1943             *val3o = val3i;
1944             break;
1945
1946         case 4 : // 1 >  2, 2 <= 3, 3 <= 1  => 1 3 2
1947             *idx1o = idx1i;
1948             *idx2o = idx3i;
1949             *idx3o = idx2i;
1950             *val1o = val1i;
1951             *val2o = val3i;
1952             *val3o = val2i;
1953             break;
1954
1955         case 2 : // 1 <= 2, 2 > 3, 3 <= 1   => 2 1 3
1956             *idx1o = idx2i;
1957             *idx2o = idx1i;
1958             *idx3o = idx3i;
1959             *val1o = val2i;
1960             *val2o = val1i;
1961             *val3o = val3i;
1962             break;
1963
1964         case 3 : // 1 <= 2, 2 > 3, 3  > 1   => 2 3 1
1965             *idx1o = idx2i;
1966             *idx2o = idx3i;
1967             *idx3o = idx1i;
1968             *val1o = val2i;
1969             *val2o = val3i;
1970             *val3o = val1i;
1971             break;
1972
1973         case 1 : // 1 <= 2, 2 <= 3, 3 > 1   => 3 2 1
1974             *idx1o = idx3i;
1975             *idx2o = idx2i;
1976             *idx3o = idx1i;
1977             *val1o = val3i;
1978             *val2o = val2i;
1979             *val3o = val1i;
1980         break;
1981
1982         case 5 : // 1 > 2, 2 <= 3, 3 >  1   => 3 1 2
1983             *idx1o = idx3i;
1984             *idx2o = idx1i;
1985             *idx3o = idx2i;
1986             *val1o = val3i;
1987             *val2o = val1i;
1988             *val3o = val2i;
1989         break;
1990
1991         default: // the rest should NEVER happen!
1992             AliError("ERROR in Sort3!!!\n");
1993         break;
1994     }
1995 //    printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
1996 }
1997
1998 void AliTRDmcmSim::Sort6To4(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
1999                                UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
2000                                UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \
2001                                UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o)
2002 {
2003   // sorting for tracklet selection
2004
2005     UShort_t idx21s, idx22s, idx23s, dummy;
2006     UShort_t val21s, val22s, val23s;
2007     UShort_t idx23as, idx23bs;
2008     UShort_t val23as, val23bs;
2009
2010     Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
2011                  idx1o, &idx21s, &idx23as,
2012                  val1o, &val21s, &val23as);
2013
2014     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2015                  idx2o, &idx22s, &idx23bs,
2016                  val2o, &val22s, &val23bs);
2017
2018     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
2019
2020     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2021                  idx3o, idx4o, &dummy,
2022                  val3o, val4o, &dummy);
2023
2024 }
2025
2026 void AliTRDmcmSim::Sort6To2Worst(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
2027                                     UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
2028                                     UShort_t *idx5o, UShort_t *idx6o)
2029 {
2030   // sorting for tracklet selection
2031
2032     UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
2033     UShort_t val21s, val22s, val23s;
2034     UShort_t idx23as, idx23bs;
2035     UShort_t val23as, val23bs;
2036
2037     Sort3(idx1i, idx2i,   idx3i, val1i, val2i, val3i,
2038                  &dummy1, &idx21s, &idx23as,
2039                  &dummy2, &val21s, &val23as);
2040
2041     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2042                  &dummy1, &idx22s, &idx23bs,
2043                  &dummy2, &val22s, &val23bs);
2044
2045     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
2046
2047     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2048                  &dummy1, &dummy2, idx6o,
2049                  &dummy3, &dummy4, &dummy5);
2050 //    printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
2051 //            idx21s,    idx23as,    idx22s,    idx23bs,    *idx5o,    *idx6o);
2052 }
2053
2054