]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDmcmSim.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding rules
[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)
dfd03fc3 654{
0c349049 655 //
dfd03fc3 656 // Produce raw data stream from this MCM and put in buf
0c349049 657 // Returns number of words filled, or negative value
658 // with -1 * number of overflowed words
659 //
dfd03fc3 660
661 UInt_t x;
dfd03fc3 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;
b0a41e80 666 Int_t nActiveADC = 0; // number of activated ADC bits in a word
dfd03fc3 667
ecf39416 668 if( !CheckInitialized() ) return 0;
669
dfd03fc3 670 if( fFeeParam->GetRAWstoreRaw() ) {
671 adc = fADCR;
672 } else {
673 adc = fADCF;
674 }
675
676 // Produce MCM header
b0a41e80 677 x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
678
dfd03fc3 679 if (nw < maxSize) {
680 buf[nw++] = x;
f793c83d 681 //printf("\nMCM header: %X ",x);
dfd03fc3 682 }
683 else {
684 of++;
685 }
686
b0a41e80 687 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
688 // n : unused , c : ADC count, m : selected ADCs
dfd03fc3 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
b0a41e80 693 x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc
694 nActiveADC++; // number of 1 in mmm....m
dfd03fc3 695 }
696 }
b0a41e80 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
dfd03fc3 700 if (nw < maxSize) {
701 buf[nw++] = x;
f793c83d 702 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
dfd03fc3 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++ ) {
b0a41e80 715 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
dfd03fc3 716 aa = !(iAdc & 1) + 2;
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;
ecf39416 721 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
722 if (nw < maxSize) {
b0a41e80 723 buf[nw++] = x;
724 //printf("%08X ",x);
ecf39416 725 }
726 else {
b0a41e80 727 of++;
ecf39416 728 }
dfd03fc3 729 }
730 }
731
732 if( of != 0 ) return -of; else return nw;
733}
734
b0a41e80 735Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
987ba9a3 736{
737 //
b0a41e80 738 // Produce tracklet data stream from this MCM and put in buf
987ba9a3 739 // Returns number of words filled, or negative value
740 // with -1 * number of overflowed words
741 //
742
987ba9a3 743 Int_t nw = 0; // Number of written words
744 Int_t of = 0; // Number of overflowed words
b0a41e80 745
987ba9a3 746 if( !CheckInitialized() ) return 0;
747
b0a41e80 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
987ba9a3 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++;
987ba9a3 756 }
b0a41e80 757
758 if( of != 0 ) return -of; else return nw;
759}
987ba9a3 760
b0a41e80 761void AliTRDmcmSim::Filter()
762{
763 //
764 // Filter the raw ADC values. The active filter stages and their
765 // parameters are taken from AliTRDtrapConfig.
766 // The raw data is stored separate from the filtered data. Thus,
767 // it is possible to run the filters on a set of raw values
768 // sequentially for parameter tuning.
769 //
987ba9a3 770
b0a41e80 771 if( !CheckInitialized() ) {
772 AliError("got called before initialization! Nothing done!");
773 return;
987ba9a3 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.
787}
987ba9a3 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).
987ba9a3 793
b0a41e80 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?
987ba9a3 797
b0a41e80 798 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
799 fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
987ba9a3 800}
801
b0a41e80 802UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
1d93b218 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.
1d93b218 807
b0a41e80 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
1d93b218 812
b0a41e80 813 UShort_t accumulatorShifted;
814 Int_t correction;
815 UShort_t inpAdd;
816
817 inpAdd = value + fpnp;
1d93b218 818
b0a41e80 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
1d93b218 827 }
828
b0a41e80 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 }
1d93b218 839}
840
b0a41e80 841void AliTRDmcmSim::FilterPedestal()
dfd03fc3 842{
0c349049 843 //
b0a41e80 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 }
b0a41e80 857}
858
859void AliTRDmcmSim::FilterGainInit()
860{
861 // Initializes the gain filter. In this case, only threshold
862 // counters are reset.
dfd03fc3 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.
23200400 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]++;
dfd03fc3 902 }
b0a41e80 903
904 return value;
dfd03fc3 905}
906
dfd03fc3 907void AliTRDmcmSim::FilterGain()
908{
b0a41e80 909 // Read data from fADCF and apply gain filter.
0c349049 910
b0a41e80 911 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
912 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
913 fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
914 }
915 }
dfd03fc3 916}
917
b0a41e80 918void AliTRDmcmSim::FilterTailInit(Int_t baseline)
dfd03fc3 919{
b0a41e80 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}
dfd03fc3 950
b0a41e80 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;
1005 }
1006}
dfd03fc3 1007
b0a41e80 1008void AliTRDmcmSim::FilterTail()
1009{
1010 // Apply tail filter
dfd03fc3 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]);
dfd03fc3 1015 }
dfd03fc3 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
ecf39416 1077void AliTRDmcmSim::DumpData( char *f, char *target )
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++;
1179}
dfd03fc3 1180
b0a41e80 1181void AliTRDmcmSim::CalcFitreg()
1182{
1183 // Preprocessing.
1184 // Detect the hits and fill the fit registers.
1185 // Requires 12-bit data from fADCF which means Filter()
1186 // has to be called before even if all filters are bypassed.
1187
1188 //???
1189 // TRAP parameters:
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;
b0a41e80 1251 }
1252 else
ab9f7002 1253 qTotal[adcch] = 0; //jkl
4ff7ed2b 1254 AliDebug(10,Form("ch %2d qTotal %5d",adcch, qTotal[adcch]));
b0a41e80 1255 }
dfd03fc3 1256
b0a41e80 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++;
1270 }
1271 adcch++;
1272 }
dfd03fc3 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--;
1286 }
dfd03fc3 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;
dfd03fc3 1293
b0a41e80 1294 found = 0;
1295 for (adcch = 0; adcch < 19; adcch++)
ab9f7002 1296 if (qTotal[adcch] > 0) found++;
b0a41e80 1297 // NOT READY
1298
1299 if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1300 {
1301 if (marked[4] == marked[5]) marked[5] = 19;
1302 for (found=0; found<6; found++)
1303 {
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]));
b0a41e80 1306 }
dfd03fc3 1307
b0a41e80 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));
b0a41e80 1321 }
1322 if (worse2 < 19)
1323 {
ab9f7002 1324 qTotal[worse2] = 0;
4ff7ed2b 1325 AliDebug(10,Form("Kill ch %d\n",worse2));
b0a41e80 1326 }
1327 }
1328
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);
b0a41e80 1408 }
dfd03fc3 1409 }
1410 }
1411}
1412
b0a41e80 1413void AliTRDmcmSim::TrackletSelection()
dfd03fc3 1414{
b0a41e80 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 }
1455 }
1456 }
1457 ntracks = 4; // cut the rest, 4 is the max
dfd03fc3 1458 }
b0a41e80 1459 // else is not necessary to sort
dfd03fc3 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;
b0a41e80 1475 }
dfd03fc3 1476 }
b0a41e80 1477 }
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}
dfd03fc3 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();
dfd03fc3 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));
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 nHits = 0;
1652 Int_t nHits0 = 0;
1653 Int_t nHits1 = 0;
40bd6ee4 1654 if (fDigitsManager) {
1655 Int_t label[30] = {0}; // up to 30 different labels possible
1656 Int_t count[30] = {0};
1657 Int_t maxIdx = -1;
1658 Int_t maxCount = 0;
1659 Int_t nLabels = 0;
1660 for (Int_t iHit = 0; iHit < fNHits; iHit++) {
1661 if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
1662 (fHits[iHit].fChannel - fFitPtr[cpu] > 1))
1663 continue;
4ff7ed2b 1664
1665 // counting contributing hits
1666 if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) &&
1667 fHits[iHit].fTimebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0))
1668 nHits0++;
1669 if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1) &&
1670 fHits[iHit].fTimebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
1671 nHits1++;
1672
40bd6ee4 1673 Int_t currLabel = fHits[iHit].fLabel;
1674 for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1675 if (currLabel == label[iLabel]) {
1676 count[iLabel]++;
1677 if (count[iLabel] > maxCount) {
1678 maxCount = count[iLabel];
1679 maxIdx = iLabel;
1680 }
1681 currLabel = 0;
1682 break;
1683 }
1684 }
1685 if (currLabel > 0) {
1686 label[nLabels++] = currLabel;
1687 }
1688 }
1689 if (maxIdx >= 0)
1690 mcLabel = label[maxIdx];
1691 }
f793c83d 1692 new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
40bd6ee4 1693 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
4ff7ed2b 1694
1695
1696 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits(fit0->fNhits + fit1->fNhits);
1697 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits0(nHits0);
1698 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1);
48e5462a 1699 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
1700 ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
b0a41e80 1701 }
dfd03fc3 1702 }
dfd03fc3 1703 }
1704}
1705
4ff7ed2b 1706Int_t AliTRDmcmSim::GetPID(Float_t q0, Float_t q1)
1707{
1708 Int_t binQ0 = (Int_t) (q0 * fgPidNBinsQ0) + 1;
1709 Int_t binQ1 = (Int_t) (q1 * fgPidNBinsQ1) + 1;
1710 binQ0 = binQ0 >= fgPidNBinsQ0 ? fgPidNBinsQ0-1 : binQ0;
1711 binQ1 = binQ1 >= fgPidNBinsQ0 ? fgPidNBinsQ0-1 : binQ1;
1712
1713 return fgPidLut[binQ0*fgPidNBinsQ1+binQ1];
1714}
1715
1716void AliTRDmcmSim::SetPIDlut(Int_t *lut, Int_t nbinsq0, Int_t nbinsq1)
1717{
1718 if (fgPidLutDelete)
1719 delete [] fgPidLut;
1720
1721 fgPidLutDelete = kFALSE;
1722 fgPidLut = lut;
1723 fgPidNBinsQ0 = nbinsq0;
1724 fgPidNBinsQ1 = nbinsq1;
1725}
1726
1727void AliTRDmcmSim::SetPIDlut(TH2F *lut)
1728{
1729 if (fgPidLutDelete)
1730 delete [] fgPidLut;
1731
1732 fgPidNBinsQ0 = lut->GetNbinsX();
1733 fgPidNBinsQ1 = lut->GetNbinsY();
1734
1735 fgPidLut = new Int_t[fgPidNBinsQ0*fgPidNBinsQ1];
1736
1737 for (Int_t ix = 0; ix < fgPidNBinsQ0; ix++) {
1738 for (Int_t iy = 0; iy < fgPidNBinsQ1; iy++) {
1739 fgPidLut[ix*fgPidNBinsQ1 + iy] = (Int_t) (256. * lut->GetBinContent(ix, iy));
1740 }
1741 }
1742
1743 fgPidLutDelete = kTRUE;
1744}
1745
1746void AliTRDmcmSim::SetPIDlutDefault()
1747{
1748 if (fgPidLutDelete )
1749 delete [] fgPidLut;
1750
1751 fgPidLutDelete = kFALSE;
1752 fgPidLut = *fgPidLutDefault;
1753 fgPidNBinsQ0 = 40;
1754 fgPidNBinsQ1 = 50;
1755}
1756
b0a41e80 1757void AliTRDmcmSim::Tracklet()
dfd03fc3 1758{
ab9f7002 1759 // Run the tracklet calculation by calling sequentially:
1760 // CalcFitreg(); TrackletSelection(); FitTracklet()
1761 // and store the tracklets
1762
b0a41e80 1763 if (!fInitialized) {
ab9f7002 1764 AliError("Called uninitialized! Nothing done!");
b0a41e80 1765 return;
dfd03fc3 1766 }
1767
b0a41e80 1768 fTrackletArray->Delete();
dfd03fc3 1769
b0a41e80 1770 CalcFitreg();
40bd6ee4 1771 if (fNHits == 0)
1772 return;
b0a41e80 1773 TrackletSelection();
1774 FitTracklet();
c8b1590d 1775}
1776
1777Bool_t AliTRDmcmSim::StoreTracklets()
1778{
40bd6ee4 1779 if (fTrackletArray->GetEntriesFast() == 0)
c8b1590d 1780 return kTRUE;
dfd03fc3 1781
b0a41e80 1782 AliRunLoader *rl = AliRunLoader::Instance();
1783 AliDataLoader *dl = 0x0;
1784 if (rl)
1785 dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1786 if (!dl) {
1787 AliError("Could not get the tracklets data loader!");
c8b1590d 1788 return kFALSE;
dfd03fc3 1789 }
b0a41e80 1790
c8b1590d 1791 TTree *trackletTree = dl->Tree();
1792 if (!trackletTree) {
1793 dl->MakeTree();
1794 trackletTree = dl->Tree();
1795 }
1796
1797 AliTRDtrackletMCM *trkl = 0x0;
1798 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1799 if (!trkbranch)
1800 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1801
1802 for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1803 trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1804 trkbranch->SetAddress(&trkl);
b0a41e80 1805// printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
c8b1590d 1806 trkbranch->Fill();
b0a41e80 1807 }
c8b1590d 1808 dl->WriteData("OVERWRITE");
1809
1810 return kTRUE;
dfd03fc3 1811}
1812
b0a41e80 1813void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
dfd03fc3 1814{
b0a41e80 1815 // write back the processed data configured by EBSF
1816 // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1817 // zero-suppressed valued are written as -1 to digits
dfd03fc3 1818
b0a41e80 1819 if (!fInitialized) {
ab9f7002 1820 AliError("Called uninitialized! Nothing done!");
b0a41e80 1821 return;
dfd03fc3 1822 }
1823
4ff7ed2b 1824// Int_t firstAdc = 0;
1825// Int_t lastAdc = fNADC - 1;
1826//
1827// while (GetCol(firstAdc) < 0)
1828// firstAdc++;
1829//
1830// while (GetCol(lastAdc) < 0)
1831// lastAdc--;
dfd03fc3 1832
4ff7ed2b 1833 Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
dfd03fc3 1834
b0a41e80 1835 if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1836 {
4ff7ed2b 1837 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
b0a41e80 1838 if (fZSM1Dim[iAdc] == 1) {
1839 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 1840 digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
b0a41e80 1841// printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin,
1842// digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1843 }
1844 }
1845 }
1846 }
1847 else {
4ff7ed2b 1848 for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
b0a41e80 1849 if (fZSM1Dim[iAdc] == 0) {
1850 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 1851 digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
b0a41e80 1852 }
1853 }
1854 else {
1855 for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
4ff7ed2b 1856 digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
b0a41e80 1857// printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1858 }
1859 }
1860 }
dfd03fc3 1861 }
b0a41e80 1862}
dfd03fc3 1863
b0a41e80 1864// help functions, to be cleaned up
1865
ab9f7002 1866UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
b0a41e80 1867{
1868 //
1869 // This function adds a and b (unsigned) and clips to
1870 // the specified number of bits.
1871 //
1872
1873 UInt_t sum = a + b;
1874 if (nbits < 32)
1875 {
1876 UInt_t maxv = (1 << nbits) - 1;;
1877 if (sum > maxv)
1878 sum = maxv;
1879 }
1880 else
1881 {
1882 if ((sum < a) || (sum < b))
1883 sum = 0xFFFFFFFF;
1884 }
1885 return sum;
dfd03fc3 1886}
1887
982869bc 1888void AliTRDmcmSim::Sort2(UShort_t idx1i, UShort_t idx2i, \
1889 UShort_t val1i, UShort_t val2i, \
1890 UShort_t *idx1o, UShort_t *idx2o, \
1891 UShort_t *val1o, UShort_t *val2o) const
dfd03fc3 1892{
ab9f7002 1893 // sorting for tracklet selection
dfd03fc3 1894
b0a41e80 1895 if (val1i > val2i)
1896 {
1897 *idx1o = idx1i;
1898 *idx2o = idx2i;
1899 *val1o = val1i;
1900 *val2o = val2i;
1901 }
1902 else
1903 {
1904 *idx1o = idx2i;
1905 *idx2o = idx1i;
1906 *val1o = val2i;
1907 *val2o = val1i;
1908 }
1909}
1910
982869bc 1911void AliTRDmcmSim::Sort3(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, \
1912 UShort_t val1i, UShort_t val2i, UShort_t val3i, \
1913 UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \
1914 UShort_t *val1o, UShort_t *val2o, UShort_t *val3o)
b0a41e80 1915{
ab9f7002 1916 // sorting for tracklet selection
1917
4ff7ed2b 1918 Int_t sel;
dfd03fc3 1919
dfd03fc3 1920
b0a41e80 1921 if (val1i > val2i) sel=4; else sel=0;
1922 if (val2i > val3i) sel=sel + 2;
1923 if (val3i > val1i) sel=sel + 1;
1924 //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1925 switch(sel)
1926 {
1927 case 6 : // 1 > 2 > 3 => 1 2 3
1928 case 0 : // 1 = 2 = 3 => 1 2 3 : in this case doesn't matter, but so is in hardware!
1929 *idx1o = idx1i;
1930 *idx2o = idx2i;
1931 *idx3o = idx3i;
1932 *val1o = val1i;
1933 *val2o = val2i;
1934 *val3o = val3i;
1935 break;
1936
1937 case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
1938 *idx1o = idx1i;
1939 *idx2o = idx3i;
1940 *idx3o = idx2i;
1941 *val1o = val1i;
1942 *val2o = val3i;
1943 *val3o = val2i;
1944 break;
1945
1946 case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
1947 *idx1o = idx2i;
1948 *idx2o = idx1i;
1949 *idx3o = idx3i;
1950 *val1o = val2i;
1951 *val2o = val1i;
1952 *val3o = val3i;
1953 break;
1954
1955 case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
1956 *idx1o = idx2i;
1957 *idx2o = idx3i;
1958 *idx3o = idx1i;
1959 *val1o = val2i;
1960 *val2o = val3i;
1961 *val3o = val1i;
1962 break;
1963
1964 case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
1965 *idx1o = idx3i;
1966 *idx2o = idx2i;
1967 *idx3o = idx1i;
1968 *val1o = val3i;
1969 *val2o = val2i;
1970 *val3o = val1i;
1971 break;
1972
1973 case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
1974 *idx1o = idx3i;
1975 *idx2o = idx1i;
1976 *idx3o = idx2i;
1977 *val1o = val3i;
1978 *val2o = val1i;
1979 *val3o = val2i;
1980 break;
1981
1982 default: // the rest should NEVER happen!
40bd6ee4 1983 AliError("ERROR in Sort3!!!\n");
b0a41e80 1984 break;
1985 }
1986// printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
1987}
dfd03fc3 1988
982869bc 1989void AliTRDmcmSim::Sort6To4(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, UShort_t idx4i, UShort_t idx5i, UShort_t idx6i, \
1990 UShort_t val1i, UShort_t val2i, UShort_t val3i, UShort_t val4i, UShort_t val5i, UShort_t val6i, \
1991 UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \
1992 UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o)
b0a41e80 1993{
ab9f7002 1994 // sorting for tracklet selection
dfd03fc3 1995
982869bc 1996 UShort_t idx21s, idx22s, idx23s, dummy;
1997 UShort_t val21s, val22s, val23s;
1998 UShort_t idx23as, idx23bs;
1999 UShort_t val23as, val23bs;
dfd03fc3 2000
b0a41e80 2001 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
2002 idx1o, &idx21s, &idx23as,
2003 val1o, &val21s, &val23as);
dfd03fc3 2004
b0a41e80 2005 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2006 idx2o, &idx22s, &idx23bs,
2007 val2o, &val22s, &val23bs);
2008
2009 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
2010
2011 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2012 idx3o, idx4o, &dummy,
2013 val3o, val4o, &dummy);
dfd03fc3 2014
dfd03fc3 2015}
2016
982869bc 2017void AliTRDmcmSim::Sort6To2Worst(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, UShort_t idx4i, UShort_t idx5i, UShort_t idx6i, \
2018 UShort_t val1i, UShort_t val2i, UShort_t val3i, UShort_t val4i, UShort_t val5i, UShort_t val6i, \
2019 UShort_t *idx5o, UShort_t *idx6o)
b0a41e80 2020{
ab9f7002 2021 // sorting for tracklet selection
1d93b218 2022
982869bc 2023 UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
2024 UShort_t val21s, val22s, val23s;
2025 UShort_t idx23as, idx23bs;
2026 UShort_t val23as, val23bs;
1d93b218 2027
b0a41e80 2028 Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
2029 &dummy1, &idx21s, &idx23as,
2030 &dummy2, &val21s, &val23as);
1d93b218 2031
b0a41e80 2032 Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2033 &dummy1, &idx22s, &idx23bs,
2034 &dummy2, &val22s, &val23bs);
1d93b218 2035
b0a41e80 2036 Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
b65e5048 2037
b0a41e80 2038 Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2039 &dummy1, &dummy2, idx6o,
2040 &dummy3, &dummy4, &dummy5);
2041// printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
2042// idx21s, idx23as, idx22s, idx23bs, *idx5o, *idx6o);
0d64b05f 2043}
f793c83d 2044
2045