1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
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. //
25 ///////////////////////////////////////////////////////////////////////////////
27 #include <fstream> // needed for raw data dump
36 #include <TClonesArray.h>
40 #include "AliRunLoader.h"
41 #include "AliLoader.h"
42 #include "AliTRDdigit.h"
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"
57 #include "TGeoGlobalMagField.h"
59 ClassImp(AliTRDmcmSim)
61 Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
63 Float_t AliTRDmcmSim::fgChargeNorm = 65000.;
64 Int_t AliTRDmcmSim::fgAddBaseline = 0;
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 }
112 Int_t (*AliTRDmcmSim::fgPidLut) = *fgPidLutDefault;
114 //_____________________________________________________________________________
115 AliTRDmcmSim::AliTRDmcmSim() : TObject()
116 ,fInitialized(kFALSE)
127 ,fTrackletArray(NULL)
136 ,fDigitsManager(NULL)
141 ,fTailAmplShort(NULL)
146 // AliTRDmcmSim default constructor
147 // By default, nothing is initialized.
148 // It is necessary to issue Init before use.
151 AliTRDmcmSim::~AliTRDmcmSim()
154 // AliTRDmcmSim destructor
158 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
159 delete [] fADCR[iadc];
160 delete [] fADCF[iadc];
161 delete [] fZSM [iadc];
170 delete [] fGainCounterA;
171 delete [] fGainCounterB;
172 delete [] fTailAmplLong;
173 delete [] fTailAmplShort;
176 fTrackletArray->Delete();
177 delete fTrackletArray;
182 void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
185 // Initialize the class with new geometry information
186 // fADC array will be reused with filled by zero
190 fFeeParam = AliTRDfeeParam::Instance();
191 fTrapConfig = AliTRDtrapConfig::Instance();
192 fSimParam = AliTRDSimParam::Instance();
193 fCommonParam = AliTRDCommonParam::Instance();
194 fCal = AliTRDcalibDB::Instance();
195 fGeo = new AliTRDgeometry();
201 fNADC = fFeeParam->GetNadcMcm();
202 fNTimeBin = fCal->GetNumberOfTimeBins();
203 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
204 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
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];
220 fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter
221 fTailAmplLong = new UShort_t[fNADC];
222 fTailAmplShort = new UShort_t[fNADC];
224 // tracklet calculation
225 fFitReg = new FitReg_t[fNADC];
226 fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets);
228 fMCMT = new UInt_t[fMaxTracklets];
231 fInitialized = kTRUE;
236 void AliTRDmcmSim::Reset()
238 // Resets the data values and internal filter registers
239 // by re-initialising them
241 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
242 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
245 fZSM [iadc][it] = 1; // Default unread = 1
247 fZSM1Dim[iadc] = 1; // Default unread = 1
248 fGainCounterA[iadc] = 0;
249 fGainCounterB[iadc] = 0;
252 for(Int_t i = 0; i < fMaxTracklets; i++) {
256 FilterPedestalInit();
258 FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
261 void AliTRDmcmSim::SetNTimebins(Int_t ntimebins)
263 fNTimeBin = ntimebins;
264 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
268 fADCR[iadc] = new Int_t[fNTimeBin];
269 fADCF[iadc] = new Int_t[fNTimeBin];
270 fZSM [iadc] = new Int_t[fNTimeBin];
274 Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm)
276 // loads the ADC data as obtained from the digitsManager for the specified MCM
281 AliError("No Runloader given");
285 AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
287 AliError("Could not get TRDLoader");
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()) {
302 if (fNTimeBin != digits->GetNtime())
303 SetNTimebins(digits->GetNtime());
305 Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm);
307 for (Int_t ch = 0; ch < fNADC; ch++) {
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);
318 for (Int_t tb = 0; tb < fNTimeBin; tb++) {
319 if (digits->GetData(padrow,padcol, tb) < 0) {
321 fADCR[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
322 fADCF[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
325 fADCR[ch][tb] = (digits->GetData(padrow, padcol, tb) + fgAddBaseline) << fgkAddDigits;
326 fADCF[ch][tb] = (digits->GetData(padrow, padcol, tb) + fgAddBaseline) << fgkAddDigits;
340 void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
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
349 // 1: pedestal 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.
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);
362 hfp->SetStats(kFALSE);
363 hfg->SetStats(kFALSE);
364 hft->SetStats(kFALSE);
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)
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);
375 valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
378 valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
380 valueg = FilterGainNextSample(1, valuep);
383 valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
384 else if (inputTail == 1)
385 valuet = FilterTailNextSample(1, valuep);
387 valuet = FilterTailNextSample(1, valueg);
389 hfp->SetBinContent(i, valuep >> 2);
390 hfg->SetBinContent(i, valueg >> 2);
391 hft->SetBinContent(i, valuet >> 2);
394 TCanvas *c = new TCanvas;
406 Bool_t AliTRDmcmSim::CheckInitialized()
409 // Check whether object is initialized
412 if( ! fInitialized ) {
413 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
418 void AliTRDmcmSim::Print(Option_t* const option) const
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.
430 printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
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);
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]);
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);
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());
469 void AliTRDmcmSim::Draw(Option_t* const option)
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:
478 // T - plot tracklets
480 TString opt = option;
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);
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);
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);
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);
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();
533 void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
536 // Store ADC data into array of raw data
539 if( !CheckInitialized() ) return;
541 if( iadc < 0 || iadc >= fNADC ) {
542 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
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;
552 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
555 // Store ADC data into array of raw data
558 if( !CheckInitialized() ) return;
560 if( iadc < 0 || iadc >= fNADC ) {
561 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
565 fADCR[iadc][it] = adc << fgkAddDigits;
566 fADCF[iadc][it] = adc << fgkAddDigits;
569 void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
571 // Set the ADC data from an AliTRDarrayADC
574 AliError("Called uninitialized! Nothing done!");
578 fDigitsManager = digitsManager;
580 if (fNTimeBin != adcArray->GetNtime())
581 SetNTimebins(adcArray->GetNtime());
583 Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
585 // Int_t firstAdc = 0;
586 // Int_t lastAdc = fNADC-1;
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);
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);
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);
612 fADCR[iAdc][iTimeBin] = (adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) + fgAddBaseline) << fgkAddDigits;
613 fADCF[iAdc][iTimeBin] = (adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) + fgAddBaseline) << fgkAddDigits;
619 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
622 // Store ADC data into array of raw data
625 if( !CheckInitialized() ) return;
627 if( iadc < 0 || iadc >= fNADC ) {
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);
637 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
640 // Return column id of the pad for the given ADC channel
643 if( !CheckInitialized() )
646 Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
647 if (col < 0 || col >= fFeeParam->GetNcol())
653 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
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
662 Int_t nw = 0; // Number of written words
663 Int_t of = 0; // Number of overflowed words
664 Int_t rawVer = fFeeParam->GetRAWversion();
666 Int_t nActiveADC = 0; // number of activated ADC bits in a word
668 if( !CheckInitialized() ) return 0;
670 if( fFeeParam->GetRAWstoreRaw() ) {
676 // Produce MCM header
677 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
681 //printf("\nMCM header: %X ",x);
687 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
688 // n : unused , c : ADC count, m : selected ADCs
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
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 );
702 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
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
712 UInt_t aa=0, a1=0, a2=0, a3=0;
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;
732 if( of != 0 ) return -of; else return nw;
735 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
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
743 Int_t nw = 0; // Number of written words
744 Int_t of = 0; // Number of overflowed words
746 if( !CheckInitialized() ) return 0;
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
751 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
753 buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
758 if( of != 0 ) return -of; else return nw;
761 void AliTRDmcmSim::Filter()
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.
771 if( !CheckInitialized() ) {
772 AliError("got called before initialization! Nothing done!");
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
782 // Non-linearity filter not implemented.
786 // Crosstalk filter not implemented.
789 void AliTRDmcmSim::FilterPedestalInit()
791 // Initializes the pedestal filter assuming that the input has
792 // been constant for a long time (compared to the time constant).
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?
798 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
799 fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
802 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
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.
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
813 UShort_t accumulatorShifted;
817 inpAdd = value + fpnp;
819 if (fpby == 0) //??? before or after update of accumulator
822 accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF; // 10 bits
823 if (timebin == 0) // the accumulator is disabled in the drift time
825 correction = (value & 0x3FF) - accumulatorShifted;
826 fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF; // 31 bits
829 if (inpAdd <= accumulatorShifted)
833 inpAdd = inpAdd - accumulatorShifted;
841 void AliTRDmcmSim::FilterPedestal()
844 // Apply pedestal filter
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.
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]);
859 void AliTRDmcmSim::FilterGainInit()
861 // Initializes the gain filter. In this case, only threshold
862 // counters are reset.
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;
872 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
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.
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;
888 tmp = (value * fgf) >> 11;
889 if (tmp > 0xFFF) tmp = 0xFFF;
892 value = AddUintClipping(tmp, fga, 12);
894 // Update threshold counters
895 // not really useful as they are cleared with every new event
896 if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
899 fGainCounterB[adc]++;
900 else if (value >= fgta)
901 fGainCounterA[adc]++;
907 void AliTRDmcmSim::FilterGain()
909 // Read data from fADCF and apply gain filter.
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]);
918 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
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.
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
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);
933 qup = (1 - lambdaL) * (1 - lambdaS);
934 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
935 Float_t kdc = qup/qdn;
941 aout = baseline - (UShort_t) kt;
942 ql = lambdaL * (1 - lambdaS) * alphaL;
943 qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
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));
951 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
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.
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
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);
966 qup = (1 - lambdaL) * (1 - lambdaS);
967 qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
968 // Float_t kdc = qup/qdn;
975 UShort_t inpVolt = value & 0xFFF; // 12 bits
977 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
981 // add the present generator outputs
982 aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
984 // calculate the difference between the input the generated signal
986 aDiff = inpVolt - aQ;
990 // the inputs to the two generators, weighted
991 alInpv = (aDiff * alphaLong) >> 11;
993 // the new values of the registers, used next time
995 tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
996 tmp = (tmp * lambdaLong) >> 11;
997 fTailAmplLong[adc] = tmp & 0xFFF;
999 tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
1000 tmp = (tmp * lambdaShort) >> 11;
1001 fTailAmplShort[adc] = tmp & 0xFFF;
1003 // the output of the filter
1008 void AliTRDmcmSim::FilterTail()
1010 // Apply tail filter
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]);
1019 void AliTRDmcmSim::ZSMapping()
1022 // Zero Suppression Mapping implemented in TRAP chip
1024 // See detail TRAP manual "Data Indication" section:
1025 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
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
1037 Int_t **adc = fADCF;
1039 if( !CheckInitialized() ) {
1040 AliError("got called uninitialized! Nothing done!");
1044 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1045 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
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
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
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)
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;
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];
1077 void AliTRDmcmSim::DumpData( const char * const f, const char * const target )
1080 // Dump data stored (for debugging).
1081 // target should contain one or multiple of the following characters
1083 // F for filtered data
1084 // Z for zero suppression map
1085 // S Raw dat astream
1086 // other characters are simply ignored
1089 UInt_t tempbuf[1024];
1091 if( !CheckInitialized() ) return;
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 );
1098 for( Int_t t=0 ; target[t] != 0 ; t++ ) {
1099 switch( target[t] ) {
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]);
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]);
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
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]);
1146 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
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
1153 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) &&
1154 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1155 fFitReg[adc].fQ0 += qtot;
1157 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1158 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1159 fFitReg[adc].fQ1 += qtot;
1161 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1162 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
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;
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;
1181 void AliTRDmcmSim::CalcFitreg()
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.
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};
1196 //??? to be clarified:
1197 UInt_t adcMask = 0xffffffff;
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;
1204 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
1205 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)
1207 timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1208 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
1209 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)
1211 timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1213 // reset the fit registers
1215 for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
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;
1227 for (timebin = timebin1; timebin < timebin2; timebin++)
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
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) );
1242 // The accumulated charge is with the pedestal!!!
1243 qtotTemp = adcLeft + adcCentral + adcRight;
1245 (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1246 (adcLeft <= adcCentral) &&
1247 (adcCentral > adcRight) )
1248 qTotal[adcch] = qtotTemp;
1253 qTotal[adcch] = 0; //jkl
1254 AliDebug(10,Form("ch %2d qTotal %5d",adcch, qTotal[adcch]));
1260 marked[4] = 19; // invalid channel
1261 marked[5] = 19; // invalid channel
1263 while ((adcch < 16) && (found < 3))
1265 if (qTotal[adcch] > 0)
1268 marked[2*found+1]=adcch;
1277 while ((adcch > 2) && (found < 3))
1279 if (qTotal[adcch] > 0)
1281 marked[2*found]=adcch;
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++)
1295 for (adcch = 0; adcch < 19; adcch++)
1296 if (qTotal[adcch] > 0) found++;
1299 if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1301 if (marked[4] == marked[5]) marked[5] = 19;
1302 for (found=0; found<6; found++)
1304 qMarked[found] = qTotal[marked[found]] >> 4;
1305 AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found]));
1308 Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1316 // Now mask the two channels with the smallest charge
1320 AliDebug(10,Form("Kill ch %d\n",worse1));
1325 AliDebug(10,Form("Kill ch %d\n",worse2));
1329 for (adcch = 0; adcch < 19; adcch++) {
1330 if (qTotal[adcch] > 0) // the channel is marked for processing
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
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)));
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;
1347 // Calculate the center of gravity
1348 // checking for adcCentral != 0 (in case of "bad" configuration)
1349 if (adcCentral == 0)
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;
1357 // label calculation
1359 if (fDigitsManager) {
1360 Int_t label[9] = { 0 }; // up to 9 different labels possible
1361 Int_t count[9] = { 0 };
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");
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));
1381 for (Int_t iPad = 0; iPad < 3; iPad++) {
1382 if (padcol[iPad] < 0)
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]) {
1389 if (count[iLabel] > maxCount) {
1390 maxCount = count[iLabel];
1397 if (currLabel > 0) {
1398 label[nLabels++] = currLabel;
1403 mcLabel = label[maxIdx];
1406 // add the hit to the fitregister
1407 AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel);
1413 void AliTRDmcmSim::TrackletSelection()
1415 // Select up to 4 tracklet candidates from the fit registers
1416 // and assign them to the CPUs.
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
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)))
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]));
1434 for (i=0; i<ntracks;i++)
1435 AliDebug(10,Form("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]));
1439 // primitive sorting according to the number of hits
1440 for (j = 0; j < (ntracks-1); j++)
1442 for (i = j+1; i < ntracks; i++)
1444 if ( (trackletCand[j][1] < trackletCand[i][1]) ||
1445 ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
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;
1457 ntracks = 4; // cut the rest, 4 is the max
1459 // else is not necessary to sort
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++)
1464 for (i = j+1; i < ntracks; i++)
1466 if (trackletCand[j][0] < trackletCand[i][0])
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;
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]));
1487 void AliTRDmcmSim::FitTracklet()
1489 // Perform the actual tracklet fit based on the fit sums
1490 // which have been filled in the fit registers.
1492 // parameters in fitred.asm (fit program)
1493 Int_t decPlaces = 5;
1496 rndAdd = (1 << (decPlaces-1)) + 1;
1497 else if (decPlaces == 1)
1499 Int_t ndriftDp = 5; // decimal places for drift time
1500 Long64_t shift = ((Long64_t) 1 << 32);
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)));
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();
1531 bz = 0.1 * fld->SolenoidField(); // kGauss -> Tesla
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);
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
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};
1558 for (Int_t cpu = 0; cpu < 4; cpu++) {
1559 if (fFitPtr[cpu] == 31)
1561 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1565 fit0 = &fFitReg[fFitPtr[cpu] ];
1566 fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1569 mult = mult << (32 + decPlaces);
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;
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;
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
1590 temp = mult * offset;
1591 offset = temp >> 32; // take the upper 32 bits
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));
1598 AliDebug(5, Form("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i", fDetector, fRobPos, fMcmPos, slope, minslope, maxslope));
1600 temp = temp * scaleD;
1601 slope = (temp >> 32);
1602 AliDebug(5, Form("slope after scaling: %i", slope));
1605 temp = temp * scaleY;
1606 offset = (temp >> 32);
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;
1613 Bool_t rejected = kFALSE;
1614 if ((slope < minslope) || (slope > maxslope))
1617 if (rejected && GetApplyCut())
1619 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1623 if (slope > 63 || slope < -64) { // wrapping in TRAP!
1624 AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
1625 fMCMT[cpu] = 0x10001000;
1629 slope = slope & 0x7F; // 7 bit
1631 if (offset > 0xfff || offset < -0xfff)
1632 AliWarning("Overflow in offset");
1633 offset = offset & 0x1FFF; // 13 bit
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)));
1640 // qTotal = (q1 / nHits) >> 1;
1641 qTotal = GetPID(q0/length/fgChargeNorm, q1/length/fgChargeNorm);
1643 AliWarning("Overflow in charge");
1644 qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1646 // assemble and store the tracklet word
1647 fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1649 // calculate MC label
1653 if (fDigitsManager) {
1654 Int_t label[30] = {0}; // up to 30 different labels possible
1655 Int_t count[30] = {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))
1664 // counting contributing hits
1665 if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) &&
1666 fHits[iHit].fTimebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0))
1668 if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1) &&
1669 fHits[iHit].fTimebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
1672 Int_t currLabel = fHits[iHit].fLabel;
1673 for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1674 if (currLabel == label[iLabel]) {
1676 if (count[iLabel] > maxCount) {
1677 maxCount = count[iLabel];
1684 if (currLabel > 0) {
1685 label[nLabels++] = currLabel;
1689 mcLabel = label[maxIdx];
1691 new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1692 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
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);
1705 Int_t AliTRDmcmSim::GetPID(Float_t q0, Float_t q1)
1707 // get PID from accumulated charges q0 and q1
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;
1714 return fgPidLut[binQ0*fgPidNBinsQ1+binQ1];
1717 void AliTRDmcmSim::SetPIDlut(Int_t *lut, Int_t nbinsq0, Int_t nbinsq1)
1719 // set a user-defined PID LUT
1724 fgPidLutDelete = kFALSE;
1726 fgPidNBinsQ0 = nbinsq0;
1727 fgPidNBinsQ1 = nbinsq1;
1730 void AliTRDmcmSim::SetPIDlut(TH2F *lut)
1732 // set a user-defined PID LUT from a 2D histogram
1737 fgPidNBinsQ0 = lut->GetNbinsX();
1738 fgPidNBinsQ1 = lut->GetNbinsY();
1740 fgPidLut = new Int_t[fgPidNBinsQ0*fgPidNBinsQ1];
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));
1748 fgPidLutDelete = kTRUE;
1751 void AliTRDmcmSim::SetPIDlutDefault()
1753 // use the default PID LUT
1755 if (fgPidLutDelete )
1758 fgPidLutDelete = kFALSE;
1759 fgPidLut = *fgPidLutDefault;
1764 void AliTRDmcmSim::Tracklet()
1766 // Run the tracklet calculation by calling sequentially:
1767 // CalcFitreg(); TrackletSelection(); FitTracklet()
1768 // and store the tracklets
1770 if (!fInitialized) {
1771 AliError("Called uninitialized! Nothing done!");
1775 fTrackletArray->Delete();
1780 TrackletSelection();
1784 Bool_t AliTRDmcmSim::StoreTracklets()
1786 // store the found tracklets via the loader
1788 if (fTrackletArray->GetEntriesFast() == 0)
1791 AliRunLoader *rl = AliRunLoader::Instance();
1792 AliDataLoader *dl = 0x0;
1794 dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1796 AliError("Could not get the tracklets data loader!");
1800 TTree *trackletTree = dl->Tree();
1801 if (!trackletTree) {
1803 trackletTree = dl->Tree();
1806 AliTRDtrackletMCM *trkl = 0x0;
1807 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1809 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
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());
1817 dl->WriteData("OVERWRITE");
1822 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
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
1828 if (!fInitialized) {
1829 AliError("Called uninitialized! Nothing done!");
1833 // Int_t firstAdc = 0;
1834 // Int_t lastAdc = fNADC - 1;
1836 // while (GetCol(firstAdc) < 0)
1839 // while (GetCol(lastAdc) < 0)
1842 Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
1844 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
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));
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);
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);
1873 // help functions, to be cleaned up
1875 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1878 // This function adds a and b (unsigned) and clips to
1879 // the specified number of bits.
1885 UInt_t maxv = (1 << nbits) - 1;;
1891 if ((sum < a) || (sum < b))
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
1902 // sorting for tracklet selection
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)
1925 // sorting for tracklet selection
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);
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!
1946 case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
1955 case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
1964 case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
1973 case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
1982 case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
1991 default: // the rest should NEVER happen!
1992 AliError("ERROR in Sort3!!!\n");
1995 // printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
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)
2003 // sorting for tracklet selection
2005 UShort_t idx21s, idx22s, idx23s, dummy;
2006 UShort_t val21s, val22s, val23s;
2007 UShort_t idx23as, idx23bs;
2008 UShort_t val23as, val23bs;
2010 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
2011 idx1o, &idx21s, &idx23as,
2012 val1o, &val21s, &val23as);
2014 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2015 idx2o, &idx22s, &idx23bs,
2016 val2o, &val22s, &val23bs);
2018 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
2020 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2021 idx3o, idx4o, &dummy,
2022 val3o, val4o, &dummy);
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)
2030 // sorting for tracklet selection
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;
2037 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
2038 &dummy1, &idx21s, &idx23as,
2039 &dummy2, &val21s, &val23as);
2041 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2042 &dummy1, &idx22s, &idx23bs,
2043 &dummy2, &val22s, &val23bs);
2045 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
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);