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