Coding rules and warnings
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
CommitLineData
dfd03fc3 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
5eba8ada 16/* $Id$ */
17
dfd03fc3 18///////////////////////////////////////////////////////////////////////////////
19// //
20// TRD MCM (Multi Chip Module) simulator //
b0a41e80 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. //
dfd03fc3 24// //
25///////////////////////////////////////////////////////////////////////////////
26
b0a41e80 27#include <fstream> // needed for raw data dump
ecf39416 28
1d93b218 29#include <TCanvas.h>
b0a41e80 30#include <TH1F.h>
31#include <TH2F.h>
32#include <TGraph.h>
33#include <TLine.h>
ecf39416 34#include <TMath.h>
b0a41e80 35#include <TRandom.h>
36#include <TClonesArray.h>
0c349049 37
dfd03fc3 38#include "AliLog.h"
b0a41e80 39#include "AliRun.h"
40#include "AliRunLoader.h"
41#include "AliLoader.h"
42#include "AliTRDdigit.h"
0c349049 43
dfd03fc3 44#include "AliTRDfeeParam.h"
b0a41e80 45#include "AliTRDtrapConfig.h"
ecf39416 46#include "AliTRDSimParam.h"
dfd03fc3 47#include "AliTRDgeometry.h"
ecf39416 48#include "AliTRDcalibDB.h"
4cc89512 49#include "AliTRDdigitsManager.h"
b65e5048 50#include "AliTRDarrayADC.h"
40bd6ee4 51#include "AliTRDarrayDictionary.h"
1d93b218 52#include "AliTRDpadPlane.h"
52c19022 53#include "AliTRDtrackletMCM.h"
b0a41e80 54#include "AliTRDmcmSim.h"
1d93b218 55
40bd6ee4 56#include "AliMagF.h"
57#include "TGeoGlobalMagField.h"
58
dfd03fc3 59ClassImp(AliTRDmcmSim)
60
40bd6ee4 61Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
62
4ff7ed2b 63Float_t AliTRDmcmSim::fgChargeNorm = 65000.;
64Int_t AliTRDmcmSim::fgAddBaseline = 0;
65
66Int_t AliTRDmcmSim::fgPidNBinsQ0 = 40;
67Int_t AliTRDmcmSim::fgPidNBinsQ1 = 50;
68Bool_t AliTRDmcmSim::fgPidLutDelete = kFALSE;
69Int_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
112Int_t (*AliTRDmcmSim::fgPidLut) = *fgPidLutDefault;
113
dfd03fc3 114//_____________________________________________________________________________
b0a41e80 115AliTRDmcmSim::AliTRDmcmSim() : TObject()
dfd03fc3 116 ,fInitialized(kFALSE)
23200400 117 ,fMaxTracklets(-1)
b0a41e80 118 ,fDetector(-1)
6e5d4cb2 119 ,fRobPos(-1)
120 ,fMcmPos(-1)
23200400 121 ,fRow (-1)
dfd03fc3 122 ,fNADC(-1)
123 ,fNTimeBin(-1)
dfd03fc3 124 ,fADCR(NULL)
125 ,fADCF(NULL)
b0a41e80 126 ,fMCMT(NULL)
127 ,fTrackletArray(NULL)
dfd03fc3 128 ,fZSM(NULL)
129 ,fZSM1Dim(NULL)
6e5d4cb2 130 ,fFeeParam(NULL)
b0a41e80 131 ,fTrapConfig(NULL)
6e5d4cb2 132 ,fSimParam(NULL)
40bd6ee4 133 ,fCommonParam(NULL)
6e5d4cb2 134 ,fCal(NULL)
135 ,fGeo(NULL)
40bd6ee4 136 ,fDigitsManager(NULL)
b0a41e80 137 ,fPedAcc(NULL)
138 ,fGainCounterA(NULL)
139 ,fGainCounterB(NULL)
140 ,fTailAmplLong(NULL)
141 ,fTailAmplShort(NULL)
142 ,fNHits(0)
143 ,fFitReg(NULL)
dfd03fc3 144{
145 //
b0a41e80 146 // AliTRDmcmSim default constructor
dfd03fc3 147 // By default, nothing is initialized.
148 // It is necessary to issue Init before use.
149}
150
dfd03fc3 151AliTRDmcmSim::~AliTRDmcmSim()
152{
153 //
154 // AliTRDmcmSim destructor
155 //
0c349049 156
b0a41e80 157 if(fInitialized) {
dfd03fc3 158 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
16e077d0 159 delete [] fADCR[iadc];
160 delete [] fADCF[iadc];
161 delete [] fZSM [iadc];
dfd03fc3 162 }
16e077d0 163 delete [] fADCR;
164 delete [] fADCF;
165 delete [] fZSM;
166 delete [] fZSM1Dim;
1d93b218 167 delete [] fMCMT;
b0a41e80 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;
1d93b218 179 }
dfd03fc3 180}
181
b0a41e80 182void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
dfd03fc3 183{
0c349049 184 //
dfd03fc3 185 // Initialize the class with new geometry information
186 // fADC array will be reused with filled by zero
0c349049 187 //
96e6312d 188
b0a41e80 189 if (!fInitialized) {
190 fFeeParam = AliTRDfeeParam::Instance();
191 fTrapConfig = AliTRDtrapConfig::Instance();
192 fSimParam = AliTRDSimParam::Instance();
40bd6ee4 193 fCommonParam = AliTRDCommonParam::Instance();
b0a41e80 194 fCal = AliTRDcalibDB::Instance();
195 fGeo = new AliTRDgeometry();
196 }
197
198 fDetector = det;
0c349049 199 fRobPos = robPos;
200 fMcmPos = mcmPos;
dfd03fc3 201 fNADC = fFeeParam->GetNadcMcm();
ecf39416 202 fNTimeBin = fCal->GetNumberOfTimeBins();
dfd03fc3 203 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
b0a41e80 204 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
23200400 205
b0a41e80 206 if (!fInitialized) {
dfd03fc3 207 fADCR = new Int_t *[fNADC];
208 fADCF = new Int_t *[fNADC];
209 fZSM = new Int_t *[fNADC];
210 fZSM1Dim = new Int_t [fNADC];
b0a41e80 211 fGainCounterA = new UInt_t[fNADC];
212 fGainCounterB = new UInt_t[fNADC];
dfd03fc3 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 }
b0a41e80 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];
dfd03fc3 229 }
230
b0a41e80 231 fInitialized = kTRUE;
232
233 Reset();
234}
235
236void AliTRDmcmSim::Reset()
237{
238 // Resets the data values and internal filter registers
239 // by re-initialising them
240
dfd03fc3 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
b0a41e80 248 fGainCounterA[iadc] = 0;
249 fGainCounterB[iadc] = 0;
dfd03fc3 250 }
ecf39416 251
1d93b218 252 for(Int_t i = 0; i < fMaxTracklets; i++) {
253 fMCMT[i] = 0;
254 }
b0a41e80 255
256 FilterPedestalInit();
257 FilterGainInit();
258 FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
259}
1d93b218 260
4ff7ed2b 261void 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
ab9f7002 274Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm)
b0a41e80 275{
276 // loads the ADC data as obtained from the digitsManager for the specified MCM
277
64e3d742 278 Init(det, rob, mcm);
b0a41e80 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
5eba8ada 291 Bool_t retval = kTRUE;
b0a41e80 292 trdLoader->LoadDigits();
40bd6ee4 293 fDigitsManager = 0x0;
b0a41e80 294 AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
295 digMgr->SetSDigits(0);
296 digMgr->CreateArrays();
297 digMgr->ReadDigits(trdLoader->TreeD());
298 AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det);
5eba8ada 299 if (digits->HasData()) {
300 digits->Expand();
301
4ff7ed2b 302 if (fNTimeBin != digits->GetNtime())
303 SetNTimebins(digits->GetNtime());
304
5eba8ada 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);
4ff7ed2b 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 }
b0a41e80 329 }
330 }
331 }
5eba8ada 332 else
333 retval = kFALSE;
4ff7ed2b 334
b0a41e80 335 delete digMgr;
4ff7ed2b 336
337 return retval;
b0a41e80 338}
339
340void 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();
dfd03fc3 404}
405
ecf39416 406Bool_t AliTRDmcmSim::CheckInitialized()
407{
0c349049 408 //
409 // Check whether object is initialized
410 //
411
ecf39416 412 if( ! fInitialized ) {
413 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
414 }
415 return fInitialized;
416}
417
ab9f7002 418void AliTRDmcmSim::Print(Option_t* const option) const
b0a41e80 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;
c8b1590d 433 if (opt.Contains("R")) {
b0a41e80 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 }
1d93b218 441 }
442
b0a41e80 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 }
1d93b218 451 }
452
b0a41e80 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",
ab9f7002 457 iHit, fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
b0a41e80 458 }
1d93b218 459 }
1d93b218 460
b0a41e80 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 }
1d93b218 466 }
b0a41e80 467}
1d93b218 468
ab9f7002 469void AliTRDmcmSim::Draw(Option_t* const option)
b0a41e80 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 }
1d93b218 494 }
b0a41e80 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 }
1d93b218 501 }
1d93b218 502 }
b0a41e80 503 hist->Draw("colz");
1d93b218 504
b0a41e80 505 if (opt.Contains("H")) {
506 TGraph *grHits = new TGraph();
507 for (Int_t iHit = 0; iHit < fNHits; iHit++) {
508 grHits->SetPoint(iHit,
ab9f7002 509 fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256.,
510 fHits[iHit].fTimebin);
b0a41e80 511 }
512 grHits->Draw("*");
513 }
1d93b218 514
b0a41e80 515 if (opt.Contains("T")) {
516 TLine *trklLines = new TLine[4];
64e3d742 517 for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
b0a41e80 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 }
1d93b218 531}
532
ab9f7002 533void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
dfd03fc3 534{
0c349049 535 //
dfd03fc3 536 // Store ADC data into array of raw data
0c349049 537 //
dfd03fc3 538
ecf39416 539 if( !CheckInitialized() ) return;
dfd03fc3 540
541 if( iadc < 0 || iadc >= fNADC ) {
b0a41e80 542 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
dfd03fc3 543 return;
544 }
545
4ff7ed2b 546 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
b0a41e80 547 fADCR[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
548 fADCF[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
dfd03fc3 549 }
550}
551
dfd03fc3 552void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
553{
0c349049 554 //
dfd03fc3 555 // Store ADC data into array of raw data
0c349049 556 //
dfd03fc3 557
ecf39416 558 if( !CheckInitialized() ) return;
dfd03fc3 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
b0a41e80 565 fADCR[iadc][it] = adc << fgkAddDigits;
566 fADCF[iadc][it] = adc << fgkAddDigits;
567}
568
40bd6ee4 569void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
b0a41e80 570{
ab9f7002 571 // Set the ADC data from an AliTRDarrayADC
572
b0a41e80 573 if (!fInitialized) {
ab9f7002 574 AliError("Called uninitialized! Nothing done!");
b0a41e80 575 return;
576 }
577
40bd6ee4 578 fDigitsManager = digitsManager;
579
4ff7ed2b 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// }
b0a41e80 603
604 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 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);
b0a41e80 610 }
611 else {
4ff7ed2b 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);
b0a41e80 614 }
615 }
616 }
dfd03fc3 617}
618
dfd03fc3 619void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
620{
0c349049 621 //
dfd03fc3 622 // Store ADC data into array of raw data
0c349049 623 //
dfd03fc3 624
ecf39416 625 if( !CheckInitialized() ) return;
dfd03fc3 626
627 if( iadc < 0 || iadc >= fNADC ) {
dfd03fc3 628 return;
629 }
630
631 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
4ff7ed2b 632 fADCR[iadc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
633 fADCF[iadc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
dfd03fc3 634 }
635}
636
dfd03fc3 637Int_t AliTRDmcmSim::GetCol( Int_t iadc )
638{
0c349049 639 //
dfd03fc3 640 // Return column id of the pad for the given ADC channel
0c349049 641 //
642
f793c83d 643 if( !CheckInitialized() )
644 return -1;
dfd03fc3 645
a6d08b7f 646 Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
647 if (col < 0 || col >= fFeeParam->GetNcol())
648 return -1;
649 else
650 return col;
dfd03fc3 651}
652
b0a41e80 653Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
987ba9a3 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;
987ba9a3 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
b0a41e80 676 // Produce MCM header
6a04e92b 677 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
b0a41e80 678
987ba9a3 679 if (nw < maxSize) {
680 buf[nw++] = x;
f793c83d 681 //printf("\nMCM header: %X ",x);
987ba9a3 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;
f793c83d 702 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
987ba9a3 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
b0a41e80 716 aa = !(iAdc & 1) + 2;
987ba9a3 717 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
b0a41e80 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;
987ba9a3 721 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
722 if (nw < maxSize) {
b0a41e80 723 buf[nw++] = x;
724 //printf("%08X ",x);
987ba9a3 725 }
726 else {
b0a41e80 727 of++;
987ba9a3 728 }
729 }
730 }
731
732 if( of != 0 ) return -of; else return nw;
733}
734
1d93b218 735Int_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
1d93b218 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
f793c83d 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++;
1d93b218 756 }
757
758 if( of != 0 ) return -of; else return nw;
759}
760
dfd03fc3 761void AliTRDmcmSim::Filter()
762{
0c349049 763 //
b0a41e80 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.
0c349049 769 //
dfd03fc3 770
b0a41e80 771 if( !CheckInitialized() ) {
772 AliError("got called before initialization! Nothing done!");
773 return;
dfd03fc3 774 }
775
b0a41e80 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.
dfd03fc3 787}
788
b0a41e80 789void 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
802UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
dfd03fc3 803{
b0a41e80 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}
23200400 840
b0a41e80 841void AliTRDmcmSim::FilterPedestal()
842{
0c349049 843 //
dfd03fc3 844 // Apply pedestal filter
0c349049 845 //
b0a41e80 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.
dfd03fc3 851
b0a41e80 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]);
dfd03fc3 855 }
856 }
857}
858
b0a41e80 859void AliTRDmcmSim::FilterGainInit()
dfd03fc3 860{
b0a41e80 861 // Initializes the gain filter. In this case, only threshold
862 // counters are reset.
0c349049 863
b0a41e80 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 }
dfd03fc3 870}
871
b0a41e80 872UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
dfd03fc3 873{
b0a41e80 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.
dfd03fc3 878
b0a41e80 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;
dfd03fc3 884
b0a41e80 885 UInt_t tmp;
dfd03fc3 886
b0a41e80 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
907void 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]);
1d93b218 914 }
b0a41e80 915 }
916}
1d93b218 917
b0a41e80 918void 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
951UShort_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
ab9f7002 975 UShort_t inpVolt = value & 0xFFF; // 12 bits
b0a41e80 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
ab9f7002 985 if (inpVolt > aQ)
986 aDiff = inpVolt - aQ;
b0a41e80 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;
dfd03fc3 1005 }
b0a41e80 1006}
dfd03fc3 1007
b0a41e80 1008void AliTRDmcmSim::FilterTail()
1009{
1010 // Apply tail filter
fabf2e09 1011
b0a41e80 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 }
dfd03fc3 1017}
1018
dfd03fc3 1019void AliTRDmcmSim::ZSMapping()
1020{
0c349049 1021 //
dfd03fc3 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
0c349049 1026 //
dfd03fc3 1027
b0a41e80 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)
f793c83d 1035 Int_t ep = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here
ecf39416 1036
b0a41e80 1037 Int_t **adc = fADCF;
dfd03fc3 1038
b0a41e80 1039 if( !CheckInitialized() ) {
ab9f7002 1040 AliError("got called uninitialized! Nothing done!");
b0a41e80 1041 return;
1042 }
1043
1044 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1045 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
dfd03fc3 1046
ecf39416 1047 // Get ADC data currently in filter buffer
b0a41e80 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
ecf39416 1051
dfd03fc3 1052 // evaluate three conditions
0c349049 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;
dfd03fc3 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++ ) {
ecf39416 1072 fZSM1Dim[iadc] &= fZSM[iadc][it];
1073 }
1074 }
1075}
1076
36dc3337 1077void AliTRDmcmSim::DumpData( const char * const f, const char * const target )
ecf39416 1078{
0c349049 1079 //
ecf39416 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
0c349049 1087 //
1088
ecf39416 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",
b0a41e80 1095 fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector),
1096 fGeo->GetSector(fDetector), fRobPos, fMcmPos );
ecf39416 1097
4ff7ed2b 1098 for( Int_t t=0 ; target[t] != 0 ; t++ ) {
ecf39416 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");
4ff7ed2b 1139 for( Int_t i = 0 ; i < s ; i++ ) {
ecf39416 1140 of << Form(" %04x %08x\n", i, tempbuf[i]);
1141 }
dfd03fc3 1142 }
1143 }
1144}
1145
b0a41e80 1146void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label)
dfd03fc3 1147{
b0a41e80 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)))
ab9f7002 1155 fFitReg[adc].fQ0 += qtot;
b0a41e80 1156
1157 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) &&
1158 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
ab9f7002 1159 fFitReg[adc].fQ1 += qtot;
b0a41e80 1160
1161 if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) &&
1162 (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1163 {
ab9f7002 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;
b0a41e80 1170 }
1171
1172 // register hits (MC info)
ab9f7002 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;
b0a41e80 1178 fNHits++;
dfd03fc3 1179}
1180
b0a41e80 1181void AliTRDmcmSim::CalcFitreg()
dfd03fc3 1182{
b0a41e80 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:
982869bc 1190 const UShort_t lutPos[128] = { // move later to some other file
b0a41e80 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:
64e3d742 1197 UInt_t adcMask = 0xffffffff;
b0a41e80 1198
ab9f7002 1199 UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
b0a41e80 1200 Short_t ypos, fromLeft, fromRight, found;
ab9f7002 1201 UShort_t qTotal[19]; // the last is dummy
1202 UShort_t marked[6], qMarked[6], worse1, worse2;
b0a41e80 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 {
ab9f7002 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;
b0a41e80 1225 }
1226
1227 for (timebin = timebin1; timebin < timebin2; timebin++)
1228 {
ab9f7002 1229 // first find the hit candidates and store the total cluster charge in qTotal array
b0a41e80 1230 // in case of not hit store 0 there.
1231 for (adcch = 0; adcch < fNADC-2; adcch++) {
ab9f7002 1232 if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
b0a41e80 1233 {
ab9f7002 1234 adcLeft = fADCF[adcch ][timebin];
1235 adcCentral = fADCF[adcch+1][timebin];
1236 adcRight = fADCF[adcch+2][timebin];
b0a41e80 1237 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1)
ab9f7002 1238 hitQual = ( (adcLeft * adcRight) <
1239 (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
b0a41e80 1240 else
ab9f7002 1241 hitQual = 1;
b0a41e80 1242 // The accumulated charge is with the pedestal!!!
ab9f7002 1243 qtotTemp = adcLeft + adcCentral + adcRight;
1244 if ( (hitQual) &&
1245 (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1246 (adcLeft <= adcCentral) &&
1247 (adcCentral > adcRight) )
1248 qTotal[adcch] = qtotTemp;
b0a41e80 1249 else
ab9f7002 1250 qTotal[adcch] = 0;
1d93b218 1251 }
b0a41e80 1252 else
ab9f7002 1253 qTotal[adcch] = 0; //jkl
4ff7ed2b 1254 AliDebug(10,Form("ch %2d qTotal %5d",adcch, qTotal[adcch]));
1d93b218 1255 }
b0a41e80 1256
1257 fromLeft = -1;
1258 adcch = 0;
1259 found = 0;
1260 marked[4] = 19; // invalid channel
1261 marked[5] = 19; // invalid channel
ab9f7002 1262 qTotal[19] = 0;
b0a41e80 1263 while ((adcch < 16) && (found < 3))
1264 {
ab9f7002 1265 if (qTotal[adcch] > 0)
b0a41e80 1266 {
1267 fromLeft = adcch;
1268 marked[2*found+1]=adcch;
1269 found++;
1d93b218 1270 }
b0a41e80 1271 adcch++;
1d93b218 1272 }
1d93b218 1273
b0a41e80 1274 fromRight = -1;
1275 adcch = 18;
1276 found = 0;
1277 while ((adcch > 2) && (found < 3))
1278 {
ab9f7002 1279 if (qTotal[adcch] > 0)
b0a41e80 1280 {
1281 marked[2*found]=adcch;
1282 found++;
1283 fromRight = adcch;
1284 }
1285 adcch--;
1d93b218 1286 }
1d93b218 1287
4ff7ed2b 1288 AliDebug(10,Form("Fromleft=%d, Fromright=%d",fromLeft, fromRight));
b0a41e80 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++)
ab9f7002 1292 qTotal[adcch] = 0;
1d93b218 1293
b0a41e80 1294 found = 0;
1295 for (adcch = 0; adcch < 19; adcch++)
ab9f7002 1296 if (qTotal[adcch] > 0) found++;
b0a41e80 1297 // NOT READY
1d93b218 1298
b0a41e80 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 {
ab9f7002 1304 qMarked[found] = qTotal[marked[found]] >> 4;
4ff7ed2b 1305 AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found]));
1d93b218 1306 }
b0a41e80 1307
1308 Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
ab9f7002 1309 qMarked[0],
1310 qMarked[3],
1311 qMarked[4],
1312 qMarked[1],
1313 qMarked[2],
1314 qMarked[5],
b0a41e80 1315 &worse1, &worse2);
1316 // Now mask the two channels with the smallest charge
1317 if (worse1 < 19)
1318 {
ab9f7002 1319 qTotal[worse1] = 0;
4ff7ed2b 1320 AliDebug(10,Form("Kill ch %d\n",worse1));
1d93b218 1321 }
b0a41e80 1322 if (worse2 < 19)
1323 {
ab9f7002 1324 qTotal[worse2] = 0;
4ff7ed2b 1325 AliDebug(10,Form("Kill ch %d\n",worse2));
1d93b218 1326 }
1d93b218 1327 }
1d93b218 1328
b0a41e80 1329 for (adcch = 0; adcch < 19; adcch++) {
ab9f7002 1330 if (qTotal[adcch] > 0) // the channel is marked for processing
b0a41e80 1331 {
ab9f7002 1332 adcLeft = fADCF[adcch ][timebin];
1333 adcCentral = fADCF[adcch+1][timebin];
1334 adcRight = fADCF[adcch+2][timebin];
b0a41e80 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
ab9f7002 1338 Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
4ff7ed2b 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)));
b0a41e80 1342
ab9f7002 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;
f793c83d 1346
b0a41e80 1347 // Calculate the center of gravity
f793c83d 1348 // checking for adcCentral != 0 (in case of "bad" configuration)
1349 if (adcCentral == 0)
1350 continue;
ab9f7002 1351 ypos = 128*(adcLeft - adcRight) / adcCentral;
b0a41e80 1352 if (ypos < 0) ypos = -ypos;
1353 // make the correction using the LUT
ab9f7002 1354 ypos = ypos + lutPos[ypos & 0x7F];
1355 if (adcLeft > adcRight) ypos = -ypos;
40bd6ee4 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) {
5eba8ada 1377 AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
40bd6ee4 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);
4ff7ed2b 1385 AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin));
40bd6ee4 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);
1d93b218 1408 }
1d93b218 1409 }
1d93b218 1410 }
b0a41e80 1411}
1d93b218 1412
b0a41e80 1413void AliTRDmcmSim::TrackletSelection()
1414{
1415 // Select up to 4 tracklet candidates from the fit registers
1416 // and assign them to the CPUs.
1417
ab9f7002 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
b0a41e80 1420
1421 ntracks = 0;
ab9f7002 1422 for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1423 if ( (fFitReg[adcIdx].fNhits
b0a41e80 1424 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
ab9f7002 1425 (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
b0a41e80 1426 >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1427 {
ab9f7002 1428 trackletCand[ntracks][0] = adcIdx;
1429 trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
4ff7ed2b 1430 AliDebug(10,Form("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]));
b0a41e80 1431 ntracks++;
1432 };
1433
4ff7ed2b 1434 for (i=0; i<ntracks;i++)
1435 AliDebug(10,Form("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]));
b0a41e80 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 {
ab9f7002 1444 if ( (trackletCand[j][1] < trackletCand[i][1]) ||
1445 ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
b0a41e80 1446 {
1447 // swap j & i
ab9f7002 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;
b0a41e80 1454 }
1d93b218 1455 }
1d93b218 1456 }
b0a41e80 1457 ntracks = 4; // cut the rest, 4 is the max
1d93b218 1458 }
b0a41e80 1459 // else is not necessary to sort
1d93b218 1460
b0a41e80 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 {
ab9f7002 1466 if (trackletCand[j][0] < trackletCand[i][0])
b0a41e80 1467 {
1468 // swap j & i
ab9f7002 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;
1d93b218 1475 }
1d93b218 1476 }
1477 }
b0a41e80 1478 for (i = 0; i < ntracks; i++) // CPUs with tracklets.
ab9f7002 1479 fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
b0a41e80 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)
4ff7ed2b 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]));
b0a41e80 1485}
1d93b218 1486
b0a41e80 1487void 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;
4ff7ed2b 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
b0a41e80 1515
1516 // should come from trapConfig (DMEM)
1517 AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
4ff7ed2b 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()));
40bd6ee4 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());
4ff7ed2b 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
b0a41e80 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();
1d93b218 1562 }
b0a41e80 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
ab9f7002 1573 nHits = fit0->fNhits + fit1->fNhits; // number of hits
1574 sumX = fit0->fSumX + fit1->fSumX;
1575 sumX2 = fit0->fSumX2 + fit1->fSumX2;
b0a41e80 1576 denom = nHits*sumX2 - sumX*sumX;
1577
1578 mult = mult / denom; // exactly like in the TRAP program
ab9f7002 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;
b0a41e80 1583
1584 slope = nHits*sumXY - sumX * sumY;
4ff7ed2b 1585 AliDebug(5, Form("slope from fitreg: %i", slope));
b0a41e80 1586 offset = sumX2*sumY - sumX * sumXY;
1587 temp = mult * slope;
1588 slope = temp >> 32; // take the upper 32 bits
4ff7ed2b 1589 slope = -slope;
b0a41e80 1590 temp = mult * offset;
1591 offset = temp >> 32; // take the upper 32 bits
1592
4ff7ed2b 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;
b0a41e80 1596 offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1d93b218 1597
4ff7ed2b 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
40bd6ee4 1613 Bool_t rejected = kFALSE;
4ff7ed2b 1614 if ((slope < minslope) || (slope > maxslope))
40bd6ee4 1615 rejected = kTRUE;
4ff7ed2b 1616
1617 if (rejected && GetApplyCut())
b0a41e80 1618 {
1619 fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1620 }
1621 else
1622 {
4ff7ed2b 1623 if (slope > 63 || slope < -64) { // wrapping in TRAP!
40bd6ee4 1624 AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
1625 fMCMT[cpu] = 0x10001000;
1626 continue;
1627 }
b0a41e80 1628
4ff7ed2b 1629 slope = slope & 0x7F; // 7 bit
1630
40bd6ee4 1631 if (offset > 0xfff || offset < -0xfff)
b0a41e80 1632 AliWarning("Overflow in offset");
1633 offset = offset & 0x1FFF; // 13 bit
1634
4ff7ed2b 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);
b0a41e80 1642 if (qTotal > 0xff)
1643 AliWarning("Overflow in charge");
1644 qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
4ff7ed2b 1645
b0a41e80 1646 // assemble and store the tracklet word
1647 fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
40bd6ee4 1648
1649 // calculate MC label
1650 Int_t mcLabel = -1;
4ff7ed2b 1651 Int_t nHits0 = 0;
1652 Int_t nHits1 = 0;
40bd6ee4 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;
4ff7ed2b 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
40bd6ee4 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 }
f793c83d 1691 new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
40bd6ee4 1692 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
4ff7ed2b 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);
48e5462a 1698 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
1699 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
b0a41e80 1700 }
1d93b218 1701 }
1d93b218 1702 }
b0a41e80 1703}
1d93b218 1704
4ff7ed2b 1705Int_t AliTRDmcmSim::GetPID(Float_t q0, Float_t q1)
1706{
36dc3337 1707 // get PID from accumulated charges q0 and q1
1708
4ff7ed2b 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
1717void AliTRDmcmSim::SetPIDlut(Int_t *lut, Int_t nbinsq0, Int_t nbinsq1)
1718{
36dc3337 1719 // set a user-defined PID LUT
1720
4ff7ed2b 1721 if (fgPidLutDelete)
1722 delete [] fgPidLut;
1723
1724 fgPidLutDelete = kFALSE;
1725 fgPidLut = lut;
1726 fgPidNBinsQ0 = nbinsq0;
1727 fgPidNBinsQ1 = nbinsq1;
1728}
1729
1730void AliTRDmcmSim::SetPIDlut(TH2F *lut)
1731{
36dc3337 1732 // set a user-defined PID LUT from a 2D histogram
1733
4ff7ed2b 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
1751void AliTRDmcmSim::SetPIDlutDefault()
1752{
36dc3337 1753 // use the default PID LUT
1754
4ff7ed2b 1755 if (fgPidLutDelete )
1756 delete [] fgPidLut;
1757
1758 fgPidLutDelete = kFALSE;
1759 fgPidLut = *fgPidLutDefault;
1760 fgPidNBinsQ0 = 40;
1761 fgPidNBinsQ1 = 50;
1762}
1763
b0a41e80 1764void AliTRDmcmSim::Tracklet()
1765{
ab9f7002 1766 // Run the tracklet calculation by calling sequentially:
1767 // CalcFitreg(); TrackletSelection(); FitTracklet()
1768 // and store the tracklets
1769
b0a41e80 1770 if (!fInitialized) {
ab9f7002 1771 AliError("Called uninitialized! Nothing done!");
b0a41e80 1772 return;
1d93b218 1773 }
96e6312d 1774
b0a41e80 1775 fTrackletArray->Delete();
96e6312d 1776
b0a41e80 1777 CalcFitreg();
40bd6ee4 1778 if (fNHits == 0)
1779 return;
b0a41e80 1780 TrackletSelection();
1781 FitTracklet();
c8b1590d 1782}
1783
1784Bool_t AliTRDmcmSim::StoreTracklets()
1785{
36dc3337 1786 // store the found tracklets via the loader
1787
40bd6ee4 1788 if (fTrackletArray->GetEntriesFast() == 0)
c8b1590d 1789 return kTRUE;
1d93b218 1790
b0a41e80 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!");
c8b1590d 1797 return kFALSE;
1d93b218 1798 }
1d93b218 1799
c8b1590d 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);
b0a41e80 1814// printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
c8b1590d 1815 trkbranch->Fill();
1d93b218 1816 }
c8b1590d 1817 dl->WriteData("OVERWRITE");
1818
1819 return kTRUE;
b0a41e80 1820}
1d93b218 1821
b0a41e80 1822void 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
1d93b218 1827
b0a41e80 1828 if (!fInitialized) {
ab9f7002 1829 AliError("Called uninitialized! Nothing done!");
b0a41e80 1830 return;
1d93b218 1831 }
1d93b218 1832
4ff7ed2b 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--;
1d93b218 1841
4ff7ed2b 1842 Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
1d93b218 1843
b0a41e80 1844 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1845 {
4ff7ed2b 1846 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
b0a41e80 1847 if (fZSM1Dim[iAdc] == 1) {
1848 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 1849 digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
b0a41e80 1850// printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin,
1851// digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1852 }
1d93b218 1853 }
b0a41e80 1854 }
1d93b218 1855 }
52c19022 1856 else {
4ff7ed2b 1857 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
b0a41e80 1858 if (fZSM1Dim[iAdc] == 0) {
1859 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 1860 digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
b0a41e80 1861 }
1862 }
1863 else {
1864 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 1865 digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
b0a41e80 1866// printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1867 }
1868 }
52c19022 1869 }
52c19022 1870 }
b0a41e80 1871}
23200400 1872
b0a41e80 1873// help functions, to be cleaned up
1d93b218 1874
ab9f7002 1875UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
b0a41e80 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;
1d93b218 1895}
b0a41e80 1896
982869bc 1897void 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
b65e5048 1901{
ab9f7002 1902 // sorting for tracklet selection
b65e5048 1903
b0a41e80 1904 if (val1i > val2i)
b65e5048 1905 {
b0a41e80 1906 *idx1o = idx1i;
1907 *idx2o = idx2i;
1908 *val1o = val1i;
1909 *val2o = val2i;
b65e5048 1910 }
b0a41e80 1911 else
b65e5048 1912 {
b0a41e80 1913 *idx1o = idx2i;
1914 *idx2o = idx1i;
1915 *val1o = val2i;
1916 *val2o = val1i;
b65e5048 1917 }
1918}
b65e5048 1919
982869bc 1920void 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)
b65e5048 1924{
ab9f7002 1925 // sorting for tracklet selection
1926
4ff7ed2b 1927 Int_t sel;
b65e5048 1928
b65e5048 1929
b0a41e80 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)
b65e5048 1935 {
b0a41e80 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!
40bd6ee4 1992 AliError("ERROR in Sort3!!!\n");
b0a41e80 1993 break;
1994 }
1995// printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
b65e5048 1996}
b0a41e80 1997
982869bc 1998void 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)
b65e5048 2002{
ab9f7002 2003 // sorting for tracklet selection
b65e5048 2004
982869bc 2005 UShort_t idx21s, idx22s, idx23s, dummy;
2006 UShort_t val21s, val22s, val23s;
2007 UShort_t idx23as, idx23bs;
2008 UShort_t val23as, val23bs;
b0a41e80 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
b65e5048 2024}
b0a41e80 2025
982869bc 2026void 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)
0d64b05f 2029{
ab9f7002 2030 // sorting for tracklet selection
0d64b05f 2031
982869bc 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;
b0a41e80 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);
b65e5048 2046
b0a41e80 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);
0d64b05f 2052}
f793c83d 2053
2054