All known overlaps removed (M. Sitta)
[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
dfd03fc3 16///////////////////////////////////////////////////////////////////////////////
17// //
18// TRD MCM (Multi Chip Module) simulator //
19// //
20///////////////////////////////////////////////////////////////////////////////
21
ecf39416 22/* $Id$ */
23
24/*
25
26 New release on 2007/08/17
27
28AliTRDmcmSim is now stably working and zero suppression function seems ok.
29From now, the default version of raw data is set to 3 in AliTRDfeeParam.
30
31The following internal parameters were abolished because it is useless and
32made trouble:
33
34 fColOfADCbeg
35 fColOfADCend
36
37GetCol member was modified accordingly.
38
39New member function DumpData was prepared for diagnostics.
40
41ZSMapping member function was debugged. It was causing crash due to
42wrong indexing in 1 dimensional numbering. Also code was shaped up better.
43
44*/
45
46/*Semi-final version of TRD raw data simulation code with zero suppression (ZS)
47similar to TRD FEE. ZS is realized by the class group:
48
49 AliTRDfeeParam
50 AliTRDmcmSim
51 AliTRDrawData
52
53AliTRDfeeParam has been modified to have more parameters like raw data
54production version and so on. AliTRDmcmSim is new class and this is the core
55of MCM (PASA+TRAP) simulator. It has still very simple function and it will be
56another project to improve this to make it closer to the reall FEE.
57AliTRDrawData has been modified to use new class AliTRDmcmSim.
dfd03fc3 58
ecf39416 59These modifications were tested on Aug. 02 HEAD version that code itself
60compiles. I'm sure there must be still bugs and we need testing by as many as
61possible persons now. Especially it seems HLT part is impacted by problems
62because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like
63fRawVersion disappeared from AliTRDrawData).
64
65In TRD definition, we have now 4 raw data versions.
66
67 0 very old offline version (by Bogdan)
68 1 test version (no zero suppression)
69 2 final version (no zero suppression)
70 3 test version (with zero suppression)
71
72The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses
73previously existing codes. If you set this to 3, AliTRDrawData changes behavior
74to use AliTRDmcmSim with ZS.
75
76Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete.
0c349049 77However it still take time because tracklet part is not yet touched.
ecf39416 78The default raw version is 2.
79
80 Ken Oyama
81*/
82
1d93b218 83// if no histo is drawn, these are obsolete
84#include <TH1.h>
85#include <TCanvas.h>
86
87// only needed if I/O of tracklets is activated
88#include <TObject.h>
89#include <TFile.h>
90#include <TTree.h>
91#include <TSystem.h>
92
ecf39416 93#include <fstream>
0c349049 94
ecf39416 95#include <TMath.h>
0c349049 96
dfd03fc3 97#include "AliLog.h"
0c349049 98
dfd03fc3 99#include "AliTRDmcmSim.h"
100#include "AliTRDfeeParam.h"
ecf39416 101#include "AliTRDSimParam.h"
dfd03fc3 102#include "AliTRDgeometry.h"
ecf39416 103#include "AliTRDcalibDB.h"
4cc89512 104#include "AliTRDdigitsManager.h"
dfd03fc3 105
1d93b218 106// additional for new tail filter and/or tracklet
107#include "AliTRDtrapAlu.h"
108#include "AliTRDpadPlane.h"
52c19022 109#include "AliTRDtrackletMCM.h"
1d93b218 110
52c19022 111#include "AliRun.h"
112#include "AliLoader.h"
1d93b218 113
dfd03fc3 114ClassImp(AliTRDmcmSim)
115
116//_____________________________________________________________________________
117AliTRDmcmSim::AliTRDmcmSim() :TObject()
118 ,fInitialized(kFALSE)
96e6312d 119 ,fNextEvent(-1)
23200400 120 ,fMaxTracklets(-1)
6e5d4cb2 121 ,fChaId(-1)
122 ,fSector(-1)
123 ,fStack(-1)
124 ,fLayer(-1)
125 ,fRobPos(-1)
126 ,fMcmPos(-1)
127 ,fNADC(-1)
128 ,fNTimeBin(-1)
23200400 129 ,fRow (-1)
6e5d4cb2 130 ,fADCR(NULL)
131 ,fADCF(NULL)
23200400 132 ,fADCT(NULL)
133 ,fPosLUT(NULL)
134 ,fMCMT(NULL)
6e5d4cb2 135 ,fZSM(NULL)
136 ,fZSM1Dim(NULL)
dfd03fc3 137 ,fFeeParam(NULL)
ecf39416 138 ,fSimParam(NULL)
139 ,fCal(NULL)
dfd03fc3 140 ,fGeo(NULL)
6e5d4cb2 141{
142 //
143 // AliTRDmcmSim default constructor
144 //
145
146 // By default, nothing is initialized.
147 // It is necessary to issue Init before use.
148}
149
150//_____________________________________________________________________________
151AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
152 :TObject(m)
1d93b218 153 ,fInitialized(kFALSE)
96e6312d 154 ,fNextEvent(-1)
23200400 155 ,fMaxTracklets(-1)
dfd03fc3 156 ,fChaId(-1)
157 ,fSector(-1)
158 ,fStack(-1)
159 ,fLayer(-1)
160 ,fRobPos(-1)
161 ,fMcmPos(-1)
162 ,fNADC(-1)
163 ,fNTimeBin(-1)
164 ,fRow(-1)
dfd03fc3 165 ,fADCR(NULL)
166 ,fADCF(NULL)
23200400 167 ,fADCT(NULL)
168 ,fPosLUT(NULL)
169 ,fMCMT(NULL)
dfd03fc3 170 ,fZSM(NULL)
171 ,fZSM1Dim(NULL)
6e5d4cb2 172 ,fFeeParam(NULL)
173 ,fSimParam(NULL)
174 ,fCal(NULL)
175 ,fGeo(NULL)
1d93b218 176
dfd03fc3 177{
178 //
6e5d4cb2 179 // AliTRDmcmSim copy constructor
dfd03fc3 180 //
181
182 // By default, nothing is initialized.
183 // It is necessary to issue Init before use.
184}
185
186//_____________________________________________________________________________
187AliTRDmcmSim::~AliTRDmcmSim()
188{
189 //
190 // AliTRDmcmSim destructor
191 //
0c349049 192
dfd03fc3 193 if( fADCR != NULL ) {
194 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
16e077d0 195 delete [] fADCR[iadc];
196 delete [] fADCF[iadc];
1d93b218 197 delete [] fADCT[iadc];
16e077d0 198 delete [] fZSM [iadc];
dfd03fc3 199 }
16e077d0 200 delete [] fADCR;
201 delete [] fADCF;
1d93b218 202 delete [] fADCT;
16e077d0 203 delete [] fZSM;
204 delete [] fZSM1Dim;
dfd03fc3 205 }
1d93b218 206
207 if(fInitialized){
208 delete [] fPosLUT;
209 delete [] fMCMT;
210 }
211
dfd03fc3 212 delete fGeo;
6e5d4cb2 213
214}
215
216//_____________________________________________________________________________
217AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
218{
219 //
220 // Assignment operator
221 //
222
223 if (this != &m) {
224 ((AliTRDmcmSim &) m).Copy(*this);
225 }
226 return *this;
227
228}
229
230//_____________________________________________________________________________
231void AliTRDmcmSim::Copy(TObject &m) const
232{
233 //
234 // Copy function
235 //
96e6312d 236 ((AliTRDmcmSim &) m).fNextEvent = 0; //new
1d93b218 237 ((AliTRDmcmSim &) m).fMaxTracklets = 0; //new
6e5d4cb2 238 ((AliTRDmcmSim &) m).fInitialized = 0;
239 ((AliTRDmcmSim &) m).fChaId = 0;
240 ((AliTRDmcmSim &) m).fSector = 0;
241 ((AliTRDmcmSim &) m).fStack = 0;
242 ((AliTRDmcmSim &) m).fLayer = 0;
243 ((AliTRDmcmSim &) m).fRobPos = 0;
244 ((AliTRDmcmSim &) m).fMcmPos = 0;
245 ((AliTRDmcmSim &) m).fNADC = 0;
246 ((AliTRDmcmSim &) m).fNTimeBin = 0;
247 ((AliTRDmcmSim &) m).fRow = 0;
248 ((AliTRDmcmSim &) m).fADCR = 0;
249 ((AliTRDmcmSim &) m).fADCF = 0;
1d93b218 250 ((AliTRDmcmSim &) m).fADCT = 0; //new
251 ((AliTRDmcmSim &) m).fPosLUT = 0; //new
252 ((AliTRDmcmSim &) m).fMCMT = 0; //new
6e5d4cb2 253 ((AliTRDmcmSim &) m).fZSM = 0;
254 ((AliTRDmcmSim &) m).fZSM1Dim = 0;
255 ((AliTRDmcmSim &) m).fFeeParam = 0;
256 ((AliTRDmcmSim &) m).fSimParam = 0;
257 ((AliTRDmcmSim &) m).fCal = 0;
258 ((AliTRDmcmSim &) m).fGeo = 0;
259
dfd03fc3 260}
261
262//_____________________________________________________________________________
23200400 263
96e6312d 264//void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
265void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent = kFALSE ) // only for readout tree (new event)
dfd03fc3 266{
0c349049 267 //
dfd03fc3 268 // Initialize the class with new geometry information
269 // fADC array will be reused with filled by zero
0c349049 270 //
96e6312d 271
272 fNextEvent = 0;
dfd03fc3 273 fFeeParam = AliTRDfeeParam::Instance();
ecf39416 274 fSimParam = AliTRDSimParam::Instance();
275 fCal = AliTRDcalibDB::Instance();
dfd03fc3 276 fGeo = new AliTRDgeometry();
0c349049 277 fChaId = chaId;
dfd03fc3 278 fSector = fGeo->GetSector( fChaId );
053767a4 279 fStack = fGeo->GetStack( fChaId );
280 fLayer = fGeo->GetLayer( fChaId );
0c349049 281 fRobPos = robPos;
282 fMcmPos = mcmPos;
dfd03fc3 283 fNADC = fFeeParam->GetNadcMcm();
ecf39416 284 fNTimeBin = fCal->GetNumberOfTimeBins();
dfd03fc3 285 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
dfd03fc3 286
1d93b218 287 fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
288
289
290
23200400 291
96e6312d 292 if (newEvent == kTRUE) {
293 fNextEvent = 1;
294 }
23200400 295
1d93b218 296
297
dfd03fc3 298 // Allocate ADC data memory if not yet done
299 if( fADCR == NULL ) {
300 fADCR = new Int_t *[fNADC];
301 fADCF = new Int_t *[fNADC];
1d93b218 302 fADCT = new Int_t *[fNADC]; //new
dfd03fc3 303 fZSM = new Int_t *[fNADC];
304 fZSM1Dim = new Int_t [fNADC];
305 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
306 fADCR[iadc] = new Int_t[fNTimeBin];
307 fADCF[iadc] = new Int_t[fNTimeBin];
1d93b218 308 fADCT[iadc] = new Int_t[fNTimeBin]; //new
dfd03fc3 309 fZSM [iadc] = new Int_t[fNTimeBin];
310 }
311 }
312
313 // Initialize ADC data
314 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
315 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
316 fADCR[iadc][it] = 0;
317 fADCF[iadc][it] = 0;
1d93b218 318 fADCT[iadc][it] = -1; //new
dfd03fc3 319 fZSM [iadc][it] = 1; // Default unread = 1
320 }
321 fZSM1Dim[iadc] = 1; // Default unread = 1
322 }
ecf39416 323
1d93b218 324 //new:
325 fPosLUT = new Int_t[128];
326 for(Int_t i = 0; i<128; i++){
327 fPosLUT[i] = 0;
328 }
329
330 fMCMT = new UInt_t[fMaxTracklets];
331 for(Int_t i = 0; i < fMaxTracklets; i++) {
332 fMCMT[i] = 0;
333 }
334
335
dfd03fc3 336 fInitialized = kTRUE;
337}
338
339//_____________________________________________________________________________
ecf39416 340Bool_t AliTRDmcmSim::CheckInitialized()
341{
0c349049 342 //
343 // Check whether object is initialized
344 //
345
ecf39416 346 if( ! fInitialized ) {
347 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
348 }
349 return fInitialized;
350}
351
352//_____________________________________________________________________________
1d93b218 353
354
355void AliTRDmcmSim::SetPosLUT() {
053767a4 356 Double_t iHi = (Double_t)fCal->GetPRFhi();
357 Double_t iLo = (Double_t)fCal->GetPRFlo();
358 Int_t nBin = fCal->GetPRFbin();
359 Int_t iOff = fLayer * nBin;
360 Int_t kNlayer = fGeo->Nlayer();
1d93b218 361
053767a4 362 Float_t *sPRFsmp = new Float_t[nBin*kNlayer];
1d93b218 363 Double_t *sPRFlayer = new Double_t[nBin];
364
365
053767a4 366 for(Int_t i = 0; i<nBin*kNlayer; i++){
1d93b218 367
368 //printf("%f\n",fCal->GetSampledPRF()[i]);
369 sPRFsmp[i] = fCal->GetSampledPRF()[i];
370
371 }
372
373 Double_t sWidth = (iHi-iLo)/((Double_t) nBin);
374 Int_t sPad = (Int_t) (1.0/sWidth);
375
376 // get the PRF for actual layer (interpolated to ibin data-points; 61 measured)
377 for(Int_t iBin = 0; iBin < nBin; iBin++){
378 sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin];
379 }
380
381 Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5); // bin-nr. for pad-position 0
382
383 Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5); // bin-nr. for pad-position 0.5
384 bin1 = bin1 + 1;
385 bin0 = bin0 + 1; //avoid negative values in aYest (start right of symmetry center)
386 while (bin0-sPad<0) {
387 bin0 = bin0 + 1;
388 }
389 while (bin1+sPad>=nBin) {
390 bin1 = bin1 - 1;
391 }
392
393 Double_t* aYest = new Double_t[bin1-bin0+1];
394
395 /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5);
396 TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5);
397 TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5);
398 TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5);
399
400 TCanvas *c1 = new TCanvas("c1","c1",800,1000);
401 hist1->Draw();
402 TCanvas *c2 = new TCanvas("c2","c2",800,1000);
403 hist2->Draw();
404 TCanvas *c3 = new TCanvas("c3","c3",800,1000);
405 hist3->Draw();
406 TCanvas *c4 = new TCanvas("c4","c4",800,1000);
407 hist4->Draw();*/
408
409 for(Int_t iBin = bin0; iBin <= bin1; iBin++){
410 aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1
411 //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo;
412 // hist1->Fill(position,aYest[iBin-bin0]);
413 }
414
415
416
417 Double_t aY[128]; // reversed function
418
419 AliTRDtrapAlu a;
420 a.Init(1,8,0,31);
421
422 for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries;
423 Double_t yest = ((Double_t)j)/256;
424
425 Int_t iBin = 0;
426 while (yest>aYest[iBin] && iBin<(bin1-bin0)) {
427 iBin = iBin+1;
428 }
429 if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) {
430 aY[j] = 0.5; // yest too big
431 //hist2->Fill(yest,aY[j]);
432
433 }
434 else {
435 Int_t bin_d = iBin + bin0 - 1;
436 Int_t bin_u = iBin + bin0;
437 Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y
438 Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y
439 Double_t yest_d = aYest[iBin-1]; // lower estimated y
440 Double_t yest_u = aYest[iBin]; // upper estimated y
441
442 aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d;
443 //hist2->Fill(yest,aY[j]);
444
445 }
446 aY[j] = aY[j] - yest;
447 //hist3->Fill(yest,aY[j]);
448 // formatting
449 a.AssignDouble(aY[j]);
450 //a.WriteWord();
451 fPosLUT[j] = a.GetValue(); // 1+8Bit value;128 entries;LUT is steered by abs(Q(i+1)-Q(i-1))/Q(i)=COG and gives the correction to COG/2
452 //hist4->Fill(yest,fPosLUT[j]);
453
454 }
455
456
457
458 delete [] sPRFsmp;
459 delete [] sPRFlayer;
460 delete [] aYest;
461
462}
463
464
465//_____________________________________________________________________________
466Int_t* AliTRDmcmSim::GetPosLUT(){
467 return fPosLUT;
468}
469
470
471
dfd03fc3 472void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
473{
0c349049 474 //
dfd03fc3 475 // Store ADC data into array of raw data
0c349049 476 //
dfd03fc3 477
ecf39416 478 if( !CheckInitialized() ) return;
dfd03fc3 479
480 if( iadc < 0 || iadc >= fNADC ) {
481 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
482 return;
483 }
484
485 for( int it = 0 ; it < fNTimeBin ; it++ ) {
486 fADCR[iadc][it] = (Int_t)(adc[it]);
487 }
488}
489
490//_____________________________________________________________________________
491void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
492{
0c349049 493 //
dfd03fc3 494 // Store ADC data into array of raw data
0c349049 495 //
dfd03fc3 496
ecf39416 497 if( !CheckInitialized() ) return;
dfd03fc3 498
499 if( iadc < 0 || iadc >= fNADC ) {
500 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
501 return;
502 }
503
504 fADCR[iadc][it] = adc;
505}
506
507//_____________________________________________________________________________
508void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
509{
0c349049 510 //
dfd03fc3 511 // Store ADC data into array of raw data
0c349049 512 //
dfd03fc3 513
ecf39416 514 if( !CheckInitialized() ) return;
dfd03fc3 515
516 if( iadc < 0 || iadc >= fNADC ) {
517 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
518 return;
519 }
520
521 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
ecf39416 522 fADCR[iadc][it] = fSimParam->GetADCbaseline();
dfd03fc3 523 }
524}
525
526//_____________________________________________________________________________
527Int_t AliTRDmcmSim::GetCol( Int_t iadc )
528{
0c349049 529 //
dfd03fc3 530 // Return column id of the pad for the given ADC channel
0c349049 531 //
532
ecf39416 533 if( !CheckInitialized() ) return -1;
dfd03fc3 534
ecf39416 535 return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
dfd03fc3 536}
537
dfd03fc3 538//_____________________________________________________________________________
539Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
540{
0c349049 541 //
dfd03fc3 542 // Produce raw data stream from this MCM and put in buf
0c349049 543 // Returns number of words filled, or negative value
544 // with -1 * number of overflowed words
545 //
dfd03fc3 546
547 UInt_t x;
548 UInt_t iEv = 0;
549 Int_t nw = 0; // Number of written words
550 Int_t of = 0; // Number of overflowed words
551 Int_t rawVer = fFeeParam->GetRAWversion();
552 Int_t **adc;
553
ecf39416 554 if( !CheckInitialized() ) return 0;
555
dfd03fc3 556 if( fFeeParam->GetRAWstoreRaw() ) {
557 adc = fADCR;
558 } else {
559 adc = fADCF;
560 }
561
562 // Produce MCM header
563 x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
564 if (nw < maxSize) {
565 buf[nw++] = x;
566 }
567 else {
568 of++;
569 }
570
571 // Produce ADC mask
572 if( rawVer >= 3 ) {
573 x = 0;
574 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
575 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
576 x = x | (1 << iAdc);
577 }
578 }
579 if (nw < maxSize) {
580 buf[nw++] = x;
581 }
582 else {
583 of++;
584 }
585 }
586
587 // Produce ADC data. 3 timebins are packed into one 32 bits word
588 // In this version, different ADC channel will NOT share the same word
589
590 UInt_t aa=0, a1=0, a2=0, a3=0;
591
592 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
ecf39416 593 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
dfd03fc3 594 aa = !(iAdc & 1) + 2;
595 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
596 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] : 0;
597 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
598 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
ecf39416 599 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
600 if (nw < maxSize) {
601 buf[nw++] = x;
602 }
603 else {
604 of++;
605 }
dfd03fc3 606 }
607 }
608
609 if( of != 0 ) return -of; else return nw;
610}
611
dfd03fc3 612//_____________________________________________________________________________
987ba9a3 613Int_t AliTRDmcmSim::ProduceRawStreamV2( UInt_t *buf, Int_t maxSize )
614{
615 //
616 // Produce raw data stream from this MCM and put in buf
617 // Returns number of words filled, or negative value
618 // with -1 * number of overflowed words
619 //
620
621 UInt_t x;
622 UInt_t iEv = 0;
623 Int_t nw = 0; // Number of written words
624 Int_t of = 0; // Number of overflowed words
625 Int_t rawVer = fFeeParam->GetRAWversion();
626 Int_t **adc;
627 Int_t nActiveADC = 0; // number of activated ADC bits in a word
628
629 if( !CheckInitialized() ) return 0;
630
631 if( fFeeParam->GetRAWstoreRaw() ) {
632 adc = fADCR;
633 } else {
634 adc = fADCF;
635 }
636
637 // Produce MCM header
638 x = (1<<31) | ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
639 if (nw < maxSize) {
640 buf[nw++] = x;
641 //printf("\nMCM header: %X ",x);
642 }
643 else {
644 of++;
645 }
646
647 // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
648 // n : unused , c : ADC count, m : selected ADCs
649 if( rawVer >= 3 ) {
650 x = 0;
651 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
652 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
653 x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc
654 nActiveADC++; // number of 1 in mmm....m
655 }
656 }
657 x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100
658 //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
659
660 if (nw < maxSize) {
661 buf[nw++] = x;
662 //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
663 }
664 else {
665 of++;
666 }
667 }
668
669 // Produce ADC data. 3 timebins are packed into one 32 bits word
670 // In this version, different ADC channel will NOT share the same word
671
672 UInt_t aa=0, a1=0, a2=0, a3=0;
673
674 for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
675 if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
676 aa = !(iAdc & 1) + 2;
677 for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
678 a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] : 0;
679 a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
680 a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
681 x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
682 if (nw < maxSize) {
683 buf[nw++] = x;
684 //printf("%08X ",x);
685 }
686 else {
687 of++;
688 }
689 }
690 }
691
692 if( of != 0 ) return -of; else return nw;
693}
694
695//_____________________________________________________________________________
1d93b218 696Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
697{
698 //
699 // Produce tracklet data stream from this MCM and put in buf
700 // Returns number of words filled, or negative value
701 // with -1 * number of overflowed words
702 //
703
704 UInt_t x;
705 Int_t nw = 0; // Number of written words
706 Int_t of = 0; // Number of overflowed words
707
708 if( !CheckInitialized() ) return 0;
709
710 // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM
711 // fMCMT is filled continuously until no more tracklet words available
712
713 Int_t wd = 0;
714 while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
715 x = fMCMT[wd];
716 if (nw < maxSize) {
717 buf[nw++] = x;
718 }
719 else {
720 of++;
721 }
722 wd++;
723 }
724
725 if( of != 0 ) return -of; else return nw;
726}
727
728
729//_____________________________________________________________________________
dfd03fc3 730void AliTRDmcmSim::Filter()
731{
0c349049 732 //
dfd03fc3 733 // Apply digital filter
0c349049 734 //
dfd03fc3 735
ecf39416 736 if( !CheckInitialized() ) return;
dfd03fc3 737
738 // Initialize filtered data array with raw data
739 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
740 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
741 fADCF[iadc][it] = fADCR[iadc][it];
742 }
743 }
744
745 // Then apply fileters one by one to filtered data array
acc49af9 746 if( fFeeParam->IsPFon() ) FilterPedestal();
747 if( fFeeParam->IsGFon() ) FilterGain();
748 if( fFeeParam->IsTFon() ) FilterTail();
dfd03fc3 749}
750
751//_____________________________________________________________________________
752void AliTRDmcmSim::FilterPedestal()
753{
23200400 754
755
0c349049 756 //
dfd03fc3 757 // Apply pedestal filter
0c349049 758 //
dfd03fc3 759
ecf39416 760 Int_t ap = fSimParam->GetADCbaseline(); // ADC instrinsic pedestal
dfd03fc3 761 Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
6e5d4cb2 762 //Int_t tc = fFeeParam->GetPFtimeConstant(); // this makes no sense yet
dfd03fc3 763
764 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
765 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
766 fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
767 }
768 }
769}
770
771//_____________________________________________________________________________
772void AliTRDmcmSim::FilterGain()
773{
0c349049 774 //
dfd03fc3 775 // Apply gain filter (not implemented)
ecf39416 776 // Later it will be implemented because gain digital filiter will
777 // increase noise level.
0c349049 778 //
779
dfd03fc3 780}
781
782//_____________________________________________________________________________
783void AliTRDmcmSim::FilterTail()
784{
0c349049 785 //
dfd03fc3 786 // Apply exponential tail filter (Bogdan's version)
0c349049 787 //
dfd03fc3 788
789 Double_t *dtarg = new Double_t[fNTimeBin];
790 Int_t *itarg = new Int_t[fNTimeBin];
791 Int_t nexp = fFeeParam->GetTFnExp();
792 Int_t tftype = fFeeParam->GetTFtype();
793
794 switch( tftype ) {
795
796 case 0: // Exponential Filter Analog Bogdan
797 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
798 FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
799 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
800 fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
801 }
802 }
803 break;
804
805 case 1: // Exponential filter digital Bogdan
806 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
807 FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
808 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
809 fADCF[iCol][iTime] = itarg[iTime];
810 }
811 }
812 break;
813
814 case 2: // Exponential filter Marian special
815 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
816 FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
817 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
818 fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
819 }
820 }
821 break;
1d93b218 822
823 //new
824 case 3: // Exponential filter using AliTRDtrapAlu class
825 for (Int_t iCol = 0; iCol < fNADC; iCol++) {
826 FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
827 for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
828 fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
829 fADCT[iCol][iTime] = itarg[iTime]; // 12bits; to be used for tracklet; tracklet will have own container;
830 }
831 }
832 break;
833
dfd03fc3 834
835 default:
836 AliError(Form("Invalid filter type %d ! \n", tftype ));
837 break;
838 }
839
840 delete dtarg;
841 delete itarg;
842}
843
844//_____________________________________________________________________________
845void AliTRDmcmSim::ZSMapping()
846{
0c349049 847 //
dfd03fc3 848 // Zero Suppression Mapping implemented in TRAP chip
849 //
850 // See detail TRAP manual "Data Indication" section:
851 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
0c349049 852 //
dfd03fc3 853
0c349049 854 Int_t eBIS = fFeeParam->GetEBsglIndThr(); // TRAP default = 0x4 (Tis=4)
855 Int_t eBIT = fFeeParam->GetEBsumIndThr(); // TRAP default = 0x28 (Tit=40)
856 Int_t eBIL = fFeeParam->GetEBindLUT(); // TRAP default = 0xf0
857 // (lookup table accept (I2,I1,I0)=(111)
858 // or (110) or (101) or (100))
859 Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
ecf39416 860 Int_t ep = AliTRDfeeParam::GetPFeffectPedestal();
861
862 if( !CheckInitialized() ) return;
dfd03fc3 863
864 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
865 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
866
ecf39416 867 // Get ADC data currently in filter buffer
0c349049 868 Int_t ap = fADCF[iadc-1][it] - ep; // previous
869 Int_t ac = fADCF[iadc ][it] - ep; // current
870 Int_t an = fADCF[iadc+1][it] - ep; // next
ecf39416 871
dfd03fc3 872 // evaluate three conditions
0c349049 873 Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection
874 Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster
875 Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak
876
877 Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table
878 Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true
879 // and d=1 means false according to TRAP manual)
880
881 fZSM[iadc][it] &= d;
882 if( eBIN == 0 ) { // turn on neighboring ADCs
883 fZSM[iadc-1][it] &= d;
884 fZSM[iadc+1][it] &= d;
dfd03fc3 885 }
0c349049 886
dfd03fc3 887 }
888 }
889
890 // do 1 dim projection
891 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
892 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
ecf39416 893 fZSM1Dim[iadc] &= fZSM[iadc][it];
894 }
895 }
0c349049 896
ecf39416 897}
898
899//_____________________________________________________________________________
900void AliTRDmcmSim::DumpData( char *f, char *target )
901{
0c349049 902 //
ecf39416 903 // Dump data stored (for debugging).
904 // target should contain one or multiple of the following characters
905 // R for raw data
906 // F for filtered data
907 // Z for zero suppression map
908 // S Raw dat astream
909 // other characters are simply ignored
0c349049 910 //
911
ecf39416 912 UInt_t tempbuf[1024];
913
914 if( !CheckInitialized() ) return;
915
916 std::ofstream of( f, std::ios::out | std::ios::app );
917 of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
918 fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
919
920 for( int t=0 ; target[t] != 0 ; t++ ) {
921 switch( target[t] ) {
922 case 'R' :
923 case 'r' :
924 of << Form("fADCR (raw ADC data)\n");
925 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
926 of << Form(" ADC %02d: ", iadc);
927 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
928 of << Form("% 4d", fADCR[iadc][it]);
929 }
930 of << Form("\n");
931 }
932 break;
933 case 'F' :
934 case 'f' :
935 of << Form("fADCF (filtered ADC data)\n");
936 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
937 of << Form(" ADC %02d: ", iadc);
938 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
939 of << Form("% 4d", fADCF[iadc][it]);
940 }
941 of << Form("\n");
942 }
943 break;
944 case 'Z' :
945 case 'z' :
946 of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
947 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
948 of << Form(" ADC %02d: ", iadc);
949 if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed
950 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
951 if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
952 }
953 of << Form("\n");
954 }
955 break;
956 case 'S' :
957 case 's' :
958 Int_t s = ProduceRawStream( tempbuf, 1024 );
959 of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
960 of << Form(" address data\n");
961 for( int i = 0 ; i < s ; i++ ) {
962 of << Form(" %04x %08x\n", i, tempbuf[i]);
963 }
dfd03fc3 964 }
965 }
966}
967
968//_____________________________________________________________________________
0c349049 969void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
970 , Int_t n, Int_t nexp)
dfd03fc3 971{
972 //
973 // Exponential filter "analog"
974 // source will not be changed
975 //
976
977 Int_t i = 0;
978 Int_t k = 0;
979 Double_t reminder[2];
980 Double_t correction;
981 Double_t result;
982 Double_t rates[2];
983 Double_t coefficients[2];
984
985 // Initialize (coefficient = alpha, rates = lambda)
986 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
987
988 Double_t r1 = (Double_t)fFeeParam->GetTFr1();
989 Double_t r2 = (Double_t)fFeeParam->GetTFr2();
990 Double_t c1 = (Double_t)fFeeParam->GetTFc1();
991 Double_t c2 = (Double_t)fFeeParam->GetTFc2();
992
993 coefficients[0] = c1;
994 coefficients[1] = c2;
995
996 Double_t dt = 0.1;
997 rates[0] = TMath::Exp(-dt/(r1));
998 rates[1] = TMath::Exp(-dt/(r2));
999
1000 // Attention: computation order is important
1001 correction = 0.0;
1002 for (k = 0; k < nexp; k++) {
1003 reminder[k] = 0.0;
1004 }
1005
1006 for (i = 0; i < n; i++) {
1007
1008 result = ((Double_t)source[i] - correction); // no rescaling
1009 target[i] = result;
1010
1011 for (k = 0; k < nexp; k++) {
1012 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1013 }
1014
1015 correction = 0.0;
1016 for (k = 0; k < nexp; k++) {
1017 correction += reminder[k];
1018 }
1019 }
1020}
1021
1022//_____________________________________________________________________________
0c349049 1023void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
1024 , Int_t nexp)
dfd03fc3 1025{
1026 //
1027 // Exponential filter "digital"
1028 // source will not be changed
1029 //
1030
ecf39416 1031 Int_t i = 0;
dfd03fc3 1032 Int_t fAlphaL = 0;
1033 Int_t fAlphaS = 0;
dfd03fc3 1034 Int_t fTailPed = 0;
dfd03fc3 1035 Int_t iAlphaL = 0;
1036 Int_t iAlphaS = 0;
dfd03fc3 1037
1038 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
1039 // initialize (coefficient = alpha, rates = lambda)
1040
1041 Double_t dt = 0.1;
1042 Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1043 Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1044 Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1045 Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1046
ecf39416 1047 Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
1048 Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
1049 Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
1050 Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
dfd03fc3 1051
1052 if (nexp == 1) {
1053 fAlphaL = (Int_t) (c1 * 2048.0);
1054 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
1055 }
1056 if (nexp == 2) {
1057 fAlphaL = (Int_t) (c1 * 2048.0);
1058 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
1059 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
1060 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
1061 }
1062
1063 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
1064 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
1065 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
1066 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
1067
1068 Int_t h1;
1069 Int_t h2;
1070 Int_t rem1;
1071 Int_t rem2;
1072 Int_t correction;
1073 Int_t result;
1074 Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
1075
1076 Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
1077 rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
1078 rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
1079
1080 // further initialization
1081 if ((rem1 + rem2) > 0x0FFF) {
1082 correction = 0x0FFF;
1083 }
1084 else {
1085 correction = (rem1 + rem2) & 0x0FFF;
1086 }
1087
1088 fTailPed = iFactor - correction;
1089
1090 for (i = 0; i < n; i++) {
1091
ecf39416 1092 result = (source[i] - correction);
1093 if (result < 0) { // Too much undershoot
dfd03fc3 1094 result = 0;
1095 }
1096
1097 target[i] = result;
1098
1099 h1 = (rem1 + ((iAlphaL * result) >> 11));
1100 if (h1 > 0x0FFF) {
1101 h1 = 0x0FFF;
1102 }
1103 else {
1104 h1 &= 0x0FFF;
1105 }
1106
1107 h2 = (rem2 + ((iAlphaS * result) >> 11));
1108 if (h2 > 0x0FFF) {
1109 h2 = 0x0FFF;
1110 }
1111 else {
1112 h2 &= 0x0FFF;
1113 }
1114
1115 rem1 = (iLambdaL * h1 ) >> 11;
1116 rem2 = (iLambdaS * h2 ) >> 11;
1117
1118 if ((rem1 + rem2) > 0x0FFF) {
1119 correction = 0x0FFF;
1120 }
1121 else {
1122 correction = (rem1 + rem2) & 0x0FFF;
1123 }
1124
1125 }
0c349049 1126
dfd03fc3 1127}
1128
1129//_____________________________________________________________________________
0c349049 1130void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
1131 , Int_t n)
dfd03fc3 1132{
1133 //
1134 // Exponential filter (M. Ivanov)
1135 // source will not be changed
1136 //
1137
1138 Int_t i = 0;
1139 Double_t sig1[100];
1140 Double_t sig2[100];
1141 Double_t sig3[100];
1142
1143 for (i = 0; i < n; i++) {
1144 sig1[i] = (Double_t)source[i];
1145 }
1146
1147 Float_t dt = 0.1;
1148 Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
1149 Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
1150
1151 FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
1152 FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
1153
1154 for (i = 0; i < n; i++) {
1155 target[i] = sig3[i];
1156 }
1157
1158}
1159
1160//______________________________________________________________________________
0c349049 1161void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
1162 , Double_t lambda, Int_t n)
dfd03fc3 1163{
1164 //
1165 // Special filter (M. Ivanov)
1166 //
1167
1168 Int_t i = 0;
1169 Double_t l = TMath::Exp(-lambda*0.5);
1170 Double_t in[1000];
1171 Double_t out[1000];
1172
1173 // Initialize in[] and out[] goes 0 ... 2*n+19
1174 for (i = 0; i < n*2+20; i++) {
1175 in[i] = out[i] = 0;
1176 }
1177
1178 // in[] goes 0, 1
1179 in[0] = ampin[0];
1180 in[1] = (ampin[0] + ampin[1]) * 0.5;
1181
1182 // Add charge to the end
1183 for (i = 0; i < 22; i++) {
1184 in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
1185 }
1186
1187 // Use arithmetic mean
1188 for (i = 1; i < n-1; i++) {
1189 in[2*i] = ampin[i]; // in[] goes 2, 3, ... , 2*n-4, 2*n-3
1190 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
1191 }
1192
1193 Double_t temp;
1194 out[2*n] = in[2*n];
1195 temp = 0;
1196 for (i = 2*n; i >= 0; i--) {
1197 out[i] = in[i] + temp;
1198 temp = l*(temp+in[i]);
1199 }
1200
1201 for (i = 0; i < n; i++){
1202 //ampout[i] = out[2*i+1]; // org
1203 ampout[i] = out[2*i];
1204 }
1205
1206}
1207
1208//______________________________________________________________________________
0c349049 1209void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
1210 , Double_t norm, Double_t lambda
1211 , Int_t n)
dfd03fc3 1212{
1213 //
1214 // Special filter (M. Ivanov)
1215 //
1216
1217 Int_t i = 0;
1218
1219 Double_t l = TMath::Exp(-lambda*0.5);
1220 Double_t k = l*(1.0 - norm*lambda*0.5);
1221 Double_t in[1000];
1222 Double_t out[1000];
1223
1224 // Initialize in[] and out[] goes 0 ... 2*n+19
1225 for (i = 0; i < n*2+20; i++) {
1226 in[i] = out[i] = 0;
1227 }
1228
1229 // in[] goes 0, 1
1230 in[0] = ampin[0];
1231 in[1] = (ampin[0]+ampin[1])*0.5;
1232
1233 // Add charge to the end
1234 for (i =-2; i < 22; i++) {
1235 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
1236 in[2*(n-1)+i] = ampin[n-1];
1237 }
1238
1239 for (i = 1; i < n-2; i++) {
1240 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
1241 in[2*i] = ampin[i];
1242 in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
1243 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
1244 }
1245
1246 Double_t temp;
1247 out[0] = in[0];
1248 temp = in[0];
1249 for (i = 1; i <= 2*n; i++) {
1250 out[i] = in[i] + (k-l)*temp;
1251 temp = in[i] + k *temp;
1252 }
1253
1254 for (i = 0; i < n; i++) {
1255 //ampout[i] = out[2*i+1]; // org
1256 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
1257 ampout[i] = TMath::Max(out[2*i],0.0);
1258 }
1259}
1260
1d93b218 1261
1262//_____________________________________________________________________________________
1263//the following filter uses AliTRDtrapAlu-class
1264
1265void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
1266 //static Int_t count = 0;
1267
1268 Double_t dt = 0.1;
1269 Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1270 Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1271 Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1272 Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1273
1274 nexp = 1;
1275
1276 //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
1277 //parameters need to be adjusted
1278 AliTRDtrapAlu lambdaL;
1279 AliTRDtrapAlu lambdaS;
1280 AliTRDtrapAlu alphaL;
1281 AliTRDtrapAlu alphaS;
1282
1283 AliTRDtrapAlu correction;
1284 AliTRDtrapAlu result;
1285 AliTRDtrapAlu bufL;
1286 AliTRDtrapAlu bufS;
1287
1288 AliTRDtrapAlu bSource;
1289
1290 lambdaL.Init(1,11);
1291 lambdaS.Init(1,11);
1292 alphaL.Init(1,11);
1293 alphaS.Init(1,11);
1294
1295 //count=count+1;
1296
1297 lambdaL.AssignDouble(TMath::Exp(-dt/r1));
1298 lambdaS.AssignDouble(TMath::Exp(-dt/r2));
1299 alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
1300 alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
1301
1302 //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
1303 correction.Init(10,2);
1304 result.Init(10,2);
1305 bufL.Init(10,2);
1306 bufS.Init(10,2);
1307 bSource.Init(10,2);
1308
1309 for(Int_t i = 0; i < n; i++) {
1310 bSource.AssignInt(source[i]);
1311 result = bSource - correction; // subtraction can produce an underflow
1312 if(result.GetSign() == kTRUE) {
1313 result.AssignInt(0);
1314 }
1315
1316 //target[i] = result.GetValuePre(); // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
1317
1318 target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2
1319
1320 //printf("target-Wert zur Zeit %d : %d",i,target[i]);
1321 //printf("\n");
1322
1323 bufL = bufL + (result * alphaL);
1324 bufL = bufL * lambdaL;
1325
1326 bufS = bufS + (result * alphaS);
1327 bufS = bufS * lambdaS; // eventually this should look like:
1328 // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
1329
1330 correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
1331 }
1332
1333
1334}
1335
1336
1337
1338
1339
1340
1341
1d93b218 1342//__________________________________________________________________________________
1343
1344
1345// in order to use the Tracklets, please first
1346// -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1347// -- set AliTRDfeeParam::fgkTFtype to 3, in order to use the new tail cancellation filter
96e6312d 1348// currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1349// -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with information about the Tracklet position (MCM, channel number)
1d93b218 1350
96e6312d 1351// The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
1d93b218 1352
1353void AliTRDmcmSim::Tracklet(){
1354 // tracklet calculation
96e6312d 1355 // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
1d93b218 1356
1357 if(!CheckInitialized()){ return; }
1358
1359 Bool_t filtered = kTRUE;
1360
1361
1362
1363 AliTRDtrapAlu data;
1364 data.Init(10,2);
1365 if(fADCT[0][0]==-1){ // check if filter was applied
1366 filtered = kFALSE;
1367 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1368 for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
1369 data.AssignInt(fADCR[iadc][iT]);
96e6312d 1370 fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
1d93b218 1371 }
1372 }
1373
1374 }
1375
96e6312d 1376 // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
1d93b218 1377 // reverse fADCT:
1378 Int_t** rev0 = new Int_t *[fNADC];
1379 Int_t** rev1 = new Int_t *[fNADC];
1380
1381 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1382 rev0[iadc] = new Int_t[fNTimeBin];
1383 rev1[iadc] = new Int_t[fNTimeBin];
1384 for( Int_t iT = 0; iT < fNTimeBin; iT++) {
1385 if( iadc <= fNADC-iadc-1 ) {
1386 rev0[iadc][iT] = fADCT[fNADC-iadc-1][iT];
1387 rev1[iadc][iT] = fADCT[iadc][iT];
1388 fADCT[iadc][iT] = rev0[iadc][iT];
1389 }
1390 else {
1391 rev0[iadc][iT] = rev1[fNADC-iadc-1][iT];
1392 fADCT[iadc][iT] = rev0[iadc][iT];
1393 }
1394 }
1395 }
1396 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1397 delete[] rev0[iadc];
1398 delete[] rev1[iadc];
1399 }
1400
1401 delete[] rev0;
1402 delete[] rev1;
1403
1404 rev0 = NULL;
1405 rev1 = NULL;
1406
1407 // get the filtered pedestal; supports only electronic tail-cancellation filter
1408 AliTRDtrapAlu filPed;
1409 Int_t ep = 0;
1410 Int_t *ieffped = new Int_t[fNTimeBin];
1411 for(Int_t iT = 0; iT < fNTimeBin; iT++){
1412 ieffped[iT] = ep;
1413 }
1414
1415 if( filtered == kTRUE ) {
1416 if( fFeeParam->IsPFon() ){
1417 ep = fFeeParam->GetPFeffectPedestal();
1418 }
1419 Int_t nexp = fFeeParam->GetTFnExp();
1420 Int_t *isource = new Int_t[fNTimeBin];
1421 filPed.Init(10,2);
1422 filPed.AssignInt(ep);
1423 Int_t epf = filPed.GetValue();
1424 for(Int_t iT = 0; iT < fNTimeBin; iT++){
1425 isource[iT] = ep;
1426 ieffped[iT] = epf;
1427 }
1428
1429 if( fFeeParam->IsTFon() ) {
1430 FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
1431 }
1432
1433 delete[] isource;
1434 }
1435
96e6312d 1436 //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
1d93b218 1437 //naming follows conventions in TRAP-manual
1438
1439
1440 Bool_t bVBY = kTRUE; // cluster-verification bypass
1441
1442 Double_t cQTParam = 0; // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
1443 AliTRDtrapAlu cQTAlu;
1444 cQTAlu.Init(1,10,0,63);
1445 cQTAlu.AssignDouble(cQTParam);
1446 Int_t cQT = cQTAlu.GetValue();
1447
1448 // linear fit
1449 Int_t tFS = fFeeParam->GetLinearFitStart(); // linear fit start
1450 Int_t tFE = fFeeParam->GetLinearFitEnd(); // linear fit stop
1451
1452 // charge accumulators
1453 Int_t tQS0 = fFeeParam->GetQacc0Start(); // start-time for charge-accumulator 0
1454 Int_t tQE0 = fFeeParam->GetQacc0End(); // stop-time for charge-accumulator 0
1455 Int_t tQS1 = fFeeParam->GetQacc1Start(); // start-time for charge-accumulator 1
1456 Int_t tQE1 = fFeeParam->GetQacc1End(); // stop-time for charge-accumulator 1
1457 // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
1458
1459 Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
1460 AliTRDtrapAlu cTHAlu;
1461 cTHAlu.Init(12,2);
1462 cTHAlu.AssignDouble(cTHParam);
1463 Int_t cTH = cTHAlu.GetValue(); // cTH used for comparison
1464
1465 struct List_t {
1466 List_t *next;
1467 Int_t iadc;
1468 Int_t value;
1469 };
1470
1471 List_t selection[7]; // list with 7 elements
1472 List_t *list = NULL;
1473 List_t *listLeft = NULL;
1474
1475 Int_t* qsum = new Int_t[fNADC];
1476
1477 // fit sums
1478 AliTRDtrapAlu qsumAlu;
1479 qsumAlu.Init(12,2); // charge sum will be 12+2 bits
1480 AliTRDtrapAlu dCOGAlu;
1481 dCOGAlu.Init(1,7,0,127); // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
1482 AliTRDtrapAlu yrawAlu;
1483 yrawAlu.Init(1,8,-1,255);
1484 AliTRDtrapAlu yAlu;
1485 yAlu.Init(1,16,-1,0xFF00); // only first 8 past-comma bits filled;additional 8 bits for accuracy;maximum 1 - 2^-8; sign is given by + or -
1486 AliTRDtrapAlu xAlu;
1487 xAlu.Init(5,8); // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
1488 AliTRDtrapAlu xxAlu;
1489 xxAlu.Init(10,0);
1490 AliTRDtrapAlu yyAlu;
1491 yyAlu.Init(1,16,0,0xFFFF); // maximum is 2^16-1; 16Bit for past-commas
1492 AliTRDtrapAlu xyAlu;
1493 xyAlu.Init(6,8);
1494 AliTRDtrapAlu XAlu;
1495 XAlu.Init(9,0);
1496 AliTRDtrapAlu XXAlu;
1497 XXAlu.Init(14,0);
1498 AliTRDtrapAlu YAlu;
1499 YAlu.Init(5,8); // 14 bit, 1 is sign-bit; therefore only 13 bit
1500 AliTRDtrapAlu YYAlu;
1501 YYAlu.Init(5,16);
1502 AliTRDtrapAlu XYAlu;
1503 XYAlu.Init(8,8); // 17 bit, 1 is sign-bit; therefore only 16 bit
1504 AliTRDtrapAlu qtruncAlu;
1505 qtruncAlu.Init(12,0);
1506 AliTRDtrapAlu QT0Alu;
1507 QT0Alu.Init(15,0);
1508 AliTRDtrapAlu QT1Alu;
1509 QT1Alu.Init(16,0);
96e6312d 1510
1511 AliTRDtrapAlu oneAlu;
1512 oneAlu.Init(1,8);
1d93b218 1513
1514
1515 AliTRDtrapAlu inverseNAlu;
96e6312d 1516 inverseNAlu.Init(1,8); // simulates the LUT for 1/N
1d93b218 1517 AliTRDtrapAlu MeanChargeAlu; // mean charge in ADC counts
1518 MeanChargeAlu.Init(8,0);
1519 AliTRDtrapAlu TotalChargeAlu;
1520 TotalChargeAlu.Init(17,8);
1521 //nr of post comma bits should be the same for inverseN and TotalCharge
1522
1523
1524 SetPosLUT(); // initialize the position correction LUT for this MCM;
1525
1526
1527 // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1528 Int_t *X = new Int_t[fNADC-2];
1529 Int_t *XX = new Int_t[fNADC-2];
1530 Int_t *Y = new Int_t[fNADC-2];
1531 Int_t *YY = new Int_t[fNADC-2];
1532 Int_t *XY = new Int_t[fNADC-2];
1533 Int_t *N = new Int_t[fNADC-2];
1534 Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1535 Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1536
1537 for (Int_t iCol = 0; iCol < fNADC-2; iCol++) {
1538
1539 // initialize fit-sums
1540 X[iCol] = 0;
1541 XX[iCol] = 0;
1542 Y[iCol] = 0;
1543 YY[iCol] = 0;
1544 XY[iCol] = 0;
1545 N[iCol] = 0;
1546 QT0[iCol] = 0;
1547 QT1[iCol] = 0;
1548 }
1549
1550
1551 filPed.Init(7,2); // convert filtered pedestal into 7+2Bits
1552
1553 for(Int_t iT = 0; iT < fNTimeBin; iT++){
1554
96e6312d 1555 if(iT<tFS || iT>=tFE) continue; // linear fit yes/no?
1d93b218 1556
1557 // reset
1558 Int_t portChannel[4] = {-1,-1,-1,-1};
1559 Int_t clusterCharge[4] = {0,0,0,0};
1560 Int_t leftCharge[4] = {0,0,0,0};
1561 Int_t centerCharge[4] = {0,0,0,0};
1562 Int_t rightCharge[4] = {0,0,0,0};
1563
1564 Int_t mark = 0;
1565
96e6312d 1566 filPed.AssignFormatted(ieffped[iT]); // no size-checking when using AssignFormatted; ieffped>=0
1d93b218 1567 filPed = filPed; // this checks the size
1568
1569 ieffped[iT] = filPed.GetValue();
1570
1571 for(Int_t i = 0; i<7; i++){
1572 selection[i].next = NULL;
1573 selection[i].iadc = -1; // value of -1: invalid adc
1574 selection[i].value = 0;
1575
1576 }
1577 // selection[0] is starting list-element; just for pointing
1578
1579 // loop over inner adc's
1580 for (Int_t iCol = 1; iCol < fNADC-1; iCol++) {
1581
1582 Int_t left = fADCT[iCol-1][iT];
1583 Int_t center = fADCT[iCol][iT];
1584 Int_t right = fADCT[iCol+1][iT];
1585
1586 Int_t sum = left + center + right; // cluster charge sum
1587 qsumAlu.AssignFormatted(sum);
1588 qsumAlu = qsumAlu; // size-checking; redundant
1589
1590 qsum[iCol] = qsumAlu.GetValue();
1591
1592 //hit detection and masking
1593 if(center>=left){
1594 if(center>right){
1595 if(qsum[iCol]>=(cTH + 3*ieffped[iT])){ // effective pedestal of all three channels must be added to cTH(+20); this is not parallel to TRAP manual; maybe cTH has to be adjusted in fFeeParam; therefore channels are not yet reduced by their pedestal
1596 mark |= 1; // marker
1597 }
1598 }
1599 }
1600 mark = mark<<1;
1601 }
1602 mark = mark>>1;
1603
1604
1605 // get selection of 6 adc's and sort,starting with greatest values
1606
1607 //read three from right side and sort (primitive sorting algorithm)
1608 Int_t i = 0; // adc number
1609 Int_t j = 1; // selection number
1610 while(i<fNADC-2 && j<=3){
1611 i = i + 1;
5a67d089 1612 if( ((mark>>(i-1)) & 1) == 1) {
1d93b218 1613 selection[j].iadc = fNADC-1-i;
1614 selection[j].value = qsum[fNADC-1-i]>>6; // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1615
1616 // insert into sorted list
1617 listLeft = &selection[0];
1618 list = listLeft->next;
1619
1620 if(list!=NULL) {
1621 while((list->next != NULL) && (selection[j].value <= list->value)){
1622 listLeft = list;
1623 list = list->next;
1624 }
1625
1626 if(selection[j].value<=list->value){
1627 selection[j].next = list->next;
1628 list->next = &selection[j];
1629 }
1630 else {
1631 listLeft->next = &selection[j];
1632 selection[j].next = list;
1633 }
1634 }
1635 else{
1636 listLeft->next = &selection[j];
1637 selection[j].next = list;
1638 }
1639
1640 j = j + 1;
1641 }
1642 }
1643
1644
1645 // read three from left side
1646 Int_t k = fNADC-2;
1647 while(k>i && j<=6) {
5a67d089 1648 if( ((mark>>(k-1)) & 1) == 1) {
1d93b218 1649 selection[j].iadc = fNADC-1-k;
1650 selection[j].value = qsum[fNADC-1-k]>>6;
1651
1652 listLeft = &selection[0];
1653 list = listLeft->next;
1654
1655 if(list!=NULL){
1656 while((list->next != NULL) && (selection[j].value <= list->value)){
1657 listLeft = list;
1658 list = list->next;
1659 }
1660
1661 if(selection[j].value<=list->value){
1662 selection[j].next = list->next;
1663 list->next = &selection[j];
1664 }
1665 else {
1666 listLeft->next = &selection[j];
1667 selection[j].next = list;
1668 }
1669 }
1670 else{
1671 listLeft->next = &selection[j];
1672 selection[j].next = list;
1673 }
1674
1675 j = j + 1;
1676 }
1677 k = k - 1;
1678 }
1679
1680 // get the four with greatest charge-sum
1681 list = &selection[0];
1682 for(i = 0; i<4; i++){
1683 if(list->next == NULL) continue;
1684 list = list->next;
1685 if(list->iadc == -1) continue;
1686 Int_t adc = list->iadc; // channel number with selected hit
1687
1688 // the following arrays contain the four chosen channels in 1 time-bin
1689 portChannel[i] = adc;
1690 clusterCharge[i] = qsum[adc];
1691 leftCharge[i] = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1692 centerCharge[i] = fADCT[adc][iT] - ieffped[iT];
1693 rightCharge[i] = fADCT[adc+1][iT] - ieffped[iT];
1694 }
1695
1696 // arithmetic unit
1697
1698 // cluster verification
1699 if(!bVBY){
1700 for(i = 0; i<4; i++){
1701 Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1702 Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1703 if (lr>=cc){
1704 portChannel[i] = -1; // set to invalid address
1705 clusterCharge[i] = 0;
1706 }
1707 }
1708 }
1709
1710 // fit-sums of valid channels
1711 // local hit position
1712 for(i = 0; i<4; i++){
1713 if (centerCharge[i] == 0) {
1714 portChannel[i] = -1;
1715 }// prevent division by 0
1716
1717 if (portChannel[i] == -1) continue;
1718
1719 Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1720
1721 Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1722 dCOG = (sign == kFALSE) ? dCOG : -dCOG; // AssignDouble doesn't allow for signed doubles
1723 dCOGAlu.AssignDouble(dCOG);
1724 Int_t iLUTpos = dCOGAlu.GetValue(); // steers position in LUT
1725
1726 dCOG = dCOG/2;
1727 yrawAlu.AssignDouble(dCOG);
1728 Int_t iCOG = yrawAlu.GetValue();
1729 Int_t y = iCOG + fPosLUT[iLUTpos % 128]; // local position in pad-units
1730 yrawAlu.AssignFormatted(y); // 0<y<1
1731 yAlu = yrawAlu; // convert to 16 past-comma bits
1732
1733 if(sign == kTRUE) yAlu.SetSign(-1); // buffer width of 9 bits; sign on real (not estimated) position
1734 xAlu.AssignInt(iT); // buffer width of 5 bits
1735
1736
1737 xxAlu = xAlu * xAlu; // buffer width of 10 bits -> fulfilled by x*x
1738
1739 yyAlu = yAlu * yAlu; // buffer width of 16 bits
1740
1741 xyAlu = xAlu * yAlu; // buffer width of 14 bits
1742
1743 Int_t adc = portChannel[i]-1; // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1744
1745 // calculate fit-sums recursively
1746 // interpretation of their bit-length is given as comment
1747
1748 // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1749
1750 XAlu.AssignFormatted(X[adc]);
1751 XAlu = XAlu + xAlu; // buffer width of 9 bits
1752 X[adc] = XAlu.GetValue();
1753
1754 XXAlu.AssignFormatted(XX[adc]);
1755 XXAlu = XXAlu + xxAlu; // buffer width of 14 bits
1756 XX[adc] = XXAlu.GetValue();
1757
1758 if (Y[adc] < 0) {
1759 YAlu.AssignFormatted(-Y[adc]); // make sure that only positive values are assigned; sign-setting must be done by hand
1760 YAlu.SetSign(-1);
1761 }
1762 else {
1763 YAlu.AssignFormatted(Y[adc]);
1764 YAlu.SetSign(1);
1765 }
1766
1767 YAlu = YAlu + yAlu; // buffer width of 14 bits (8 past-comma);
1768 Y[adc] = YAlu.GetSignedValue();
1769
1770 YYAlu.AssignFormatted(YY[adc]);
1771 YYAlu = YYAlu + yyAlu; // buffer width of 21 bits (16 past-comma)
1772 YY[adc] = YYAlu.GetValue();
1773
1774 if (XY[adc] < 0) {
1775 XYAlu.AssignFormatted(-XY[adc]);
1776 XYAlu.SetSign(-1);
1777 }
1778 else {
1779 XYAlu.AssignFormatted(XY[adc]);
1780 XYAlu.SetSign(1);
1781 }
1782
1783 XYAlu = XYAlu + xyAlu; // buffer allows 17 bits (8 past-comma)
1784 XY[adc] = XYAlu.GetSignedValue();
1785
1786 N[adc] = N[adc] + 1;
1787
1788
1789 // accumulated charge
1790 qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1791 qtruncAlu = qsumAlu;
1792
1793 if(iT>=tQS0 && iT<=tQE0){
1794 QT0Alu.AssignFormatted(QT0[adc]);
1795 QT0Alu = QT0Alu + qtruncAlu;
1796 QT0[adc] = QT0Alu.GetValue();
1797 //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1798 }
1799
1800 if(iT>=tQS1 && iT<=tQE1){
1801 QT1Alu.AssignFormatted(QT1[adc]);
1802 QT1Alu = QT1Alu + qtruncAlu;
1803 QT1[adc] = QT1Alu.GetValue();
1804 //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1805 }
1806 }// i
1807
1808 // remapping is done!!
1809
1810 }//iT
1811
1812
1813
1814 // tracklet-assembly
1815
1816 // put into AliTRDfeeParam and take care that values are in proper range
1817 const Int_t cTCL = 1; // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??)
96e6312d 1818 const Int_t cTCT = 8; // joint number of hits; 8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
1d93b218 1819
1820 Int_t mPair = 0; // marker for possible tracklet pairs
1821 Int_t* hitSum = new Int_t[fNADC-3];
1822 // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18;
1823
1824 // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1825 for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1826 hitSum[iCol] = N[iCol] + N[iCol+1];
1827 if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1828 mPair |= 1; // mark as possible channel-pair
1829
1830 }
1831 mPair = mPair<<1;
1832 }
1833 mPair = mPair>>1;
1834
1835 List_t* selectPair = new List_t[fNADC-2]; // list with 18 elements (0..18) containing the left channel-nr and hit sums
1836 // selectPair[18] is starting list-element just for pointing
1837 for(Int_t k = 0; k<fNADC-2; k++){
1838 selectPair[k].next = NULL;
1839 selectPair[k].iadc = -1; // invalid adc
1840 selectPair[k].value = 0;
1841
1842 }
1843
1844 list = NULL;
1845 listLeft = NULL;
1846
1847 // read marker and sort according to hit-sum
1848
1849 Int_t adcL = 0; // left adc-channel-number (remapped)
1850 Int_t selNr = 0; // current number in list
1851
1852 // insert marked channels into list and sort according to hit-sum
1853 while(adcL < fNADC-3 && selNr < fNADC-3){
1854
5a67d089 1855 if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
1d93b218 1856 selectPair[selNr].iadc = adcL;
1857 selectPair[selNr].value = hitSum[adcL];
1858
1859 listLeft = &selectPair[fNADC-3];
1860 list = listLeft->next;
1861
1862 if(list!=NULL) {
1863 while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1864 listLeft = list;
1865 list = list->next;
1866 }
1867
1868 if(selectPair[selNr].value <= list->value){
1869 selectPair[selNr].next = list->next;
1870 list->next = &selectPair[selNr];
1871 }
1872 else {
1873 listLeft->next = &selectPair[selNr];
1874 selectPair[selNr].next = list;
1875 }
1876
1877 }
1878 else{
1879 listLeft->next = &selectPair[selNr];
1880 selectPair[selNr].next = list;
1881 }
1882
1883 selNr = selNr + 1;
1884 }
1885 adcL = adcL + 1;
1886 }
1887
1888 //select up to 4 channels with maximum number of hits
1889 Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1890 Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1891 list = &selectPair[fNADC-3];
1892
1893 for (Int_t i = 0; i<4; i++) {
1894 if(list->next == NULL) continue;
1895 list = list->next;
1896 if(list->iadc == -1) continue;
1897 lpairChannel[i] = list->iadc; // channel number with selected hit
1898 rpairChannel[i] = lpairChannel[i]+1;
1899 }
1900
1901 // avoid submission of double tracklets
1902 for (Int_t i = 3; i>0; i--) {
1903 for (Int_t j = i-1; j>-1; j--) {
1904 if(lpairChannel[i] == rpairChannel[j]) {
1905 lpairChannel[i] = -1;
1906 rpairChannel[i] = -1;
1907 break;
1908 }
96e6312d 1909 /* if(rpairChannel[i] == lpairChannel[j]) {
1d93b218 1910 lpairChannel[i] = -1;
1911 rpairChannel[i] = -1;
1912 break;
96e6312d 1913 }*/
1d93b218 1914 }
1915 }
1916
1917 // merging of the fit-sums of the remainig channels
1918 // assume same data-word-width as for fit-sums for 1 channel
1919 // relative scales!
1920 Int_t mADC[4];
1921 Int_t mN[4];
1922 Int_t mQT0[4];
1923 Int_t mQT1[4];
1924 Int_t mX[4];
1925 Int_t mXX[4];
1926 Int_t mY[4];
1927 Int_t mYY[4];
1928 Int_t mXY[4];
1929 Int_t mOffset[4];
1930 Int_t mSlope[4];
1931 Int_t mMeanCharge[4];
1932 Int_t inverseN = 0;
1933 Double_t invN = 0;
96e6312d 1934 Int_t one = 0;
1935
1d93b218 1936 for (Int_t i = 0; i<4; i++){
1937 mADC[i] = -1; // set to invalid number
1938 mN[i] = 0;
1939 mQT0[i] = 0;
1940 mQT1[i] = 0;
1941 mX[i] = 0;
1942 mXX[i] = 0;
1943 mY[i] = 0;
1944 mYY[i] = 0;
1945 mXY[i] = 0;
1946 mOffset[i] = 0;
1947 mSlope[i] = 0;
1948 mMeanCharge[i] = 0;
1949 }
1950
96e6312d 1951 oneAlu.AssignInt(1);
1952 one = oneAlu.GetValue(); // one with 8 past comma bits
1d93b218 1953
1954 for (Int_t i = 0; i<4; i++){
96e6312d 1955
1956
1d93b218 1957 mADC[i] = lpairChannel[i]; // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
96e6312d 1958 // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1d93b218 1959 Int_t madc = mADC[i];
1960 if (madc == -1) continue;
96e6312d 1961
1962 YAlu.AssignInt(N[rpairChannel[i]]);
1963 Int_t wpad = YAlu.GetValue(); // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
1964
1d93b218 1965 mN[i] = hitSum[madc];
1966
1967 // don't merge fit sums in case of a stand-alone tracklet (consisting of only 1 channel); in that case only left channel makes up the fit sums
1968 if (N[madc+1] == 0) {
1969 mQT0[i] = QT0[madc];
1970 mQT1[i] = QT1[madc];
1971
1972 }
1973 else {
1974
1975 // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1976
1977 mQT0[i] = QT0[madc] + QT0[madc+1];
1978 QT0Alu.AssignFormatted(mQT0[i]);
1979 QT0Alu = QT0Alu; // size-check
1980 mQT0[i] = QT0Alu.GetValue(); // write back
1981
1982 mQT1[i] = QT1[madc] + QT1[madc+1];
1983 QT1Alu.AssignFormatted(mQT1[i]);
1984 QT1Alu = QT1Alu;
1985 mQT1[i] = QT1Alu.GetValue();
1986 }
1987
1988 // calculate the mean charge in adc values; later to be replaced by electron likelihood
1989 mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1990 mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1991 // simulate LUT for 1/N; LUT is fed with the double-accurate pre-calculated value of 1/N; accuracy of entries has to be adjusted to real TRAP
1992 invN = 1.0/(mN[i]);
1993 inverseNAlu.AssignDouble(invN);
1994 inverseN = inverseNAlu.GetValue();
1995 mMeanCharge[i] = mMeanCharge[i] * inverseN; // now to be interpreted with 8 past-comma bits
1996 TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1997 TotalChargeAlu = TotalChargeAlu;
1998 MeanChargeAlu = TotalChargeAlu;
1999 mMeanCharge[i] = MeanChargeAlu.GetValue();
2000
96e6312d 2001 // this check is not necessary; it is just for efficiency reasons
1d93b218 2002 if (N[madc+1] == 0) {
2003 mX[i] = X[madc];
2004 mXX[i] = XX[madc];
2005 mY[i] = Y[madc];
2006 mXY[i] = XY[madc];
2007 mYY[i] = YY[madc];
2008 }
2009 else {
2010
2011 mX[i] = X[madc] + X[madc+1];
2012 XAlu.AssignFormatted(mX[i]);
96e6312d 2013 XAlu = XAlu;
1d93b218 2014 mX[i] = XAlu.GetValue();
2015
2016 mXX[i] = XX[madc] + XX[madc+1];
2017 XXAlu.AssignFormatted(mXX[i]);
96e6312d 2018 XXAlu = XXAlu;
1d93b218 2019 mXX[i] = XXAlu.GetValue();
2020
2021
2022 mY[i] = Y[madc] + Y[madc+1] + wpad;
2023 if (mY[i] < 0) {
2024 YAlu.AssignFormatted(-mY[i]);
2025 YAlu.SetSign(-1);
2026 }
2027 else {
2028 YAlu.AssignFormatted(mY[i]);
2029 YAlu.SetSign(1);
2030 }
2031 YAlu = YAlu;
2032 mY[i] = YAlu.GetSignedValue();
2033
96e6312d 2034 mXY[i] = XY[madc] + XY[madc+1] + X[madc+1]*one; // multiplication by one to maintain the data format
1d93b218 2035
2036 if (mXY[i] < 0) {
2037 XYAlu.AssignFormatted(-mXY[i]);
2038 XYAlu.SetSign(-1);
2039 }
2040 else {
2041 XYAlu.AssignFormatted(mXY[i]);
2042 XYAlu.SetSign(1);
2043 }
2044 XYAlu = XYAlu;
2045 mXY[i] = XYAlu.GetSignedValue();
2046
96e6312d 2047 mYY[i] = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
1d93b218 2048 if (mYY[i] < 0) {
2049 YYAlu.AssignFormatted(-mYY[i]);
2050 YYAlu.SetSign(-1);
2051 }
2052 else {
2053 YYAlu.AssignFormatted(mYY[i]);
2054 YYAlu.SetSign(1);
2055 }
2056
2057 YYAlu = YYAlu;
2058 mYY[i] = YYAlu.GetSignedValue();
2059 }
2060
2061 }
2062
2063 // calculation of offset and slope from the merged fit-sums;
2064 // YY is needed for some error measure only; still to be done
2065 // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
2066
23200400 2067 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2068 // !!important note: the offset is calculated from hits in the time bin range between tFS and tFE; it corresponds to the value at the height of the time bin tFS which does NOT need to correspond to the upper side of the drift !!
2069 // !!volume (cathode wire plane). The offset cannot be rescaled as long as it is unknown which is the first time bin that contains hits from the drift region and thus to which distance from the cathode plane tFS corresponds. !!
2070 // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below). !!
2071 // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used. !!
2072 // !!The corrections to the offset (e.g. no ExB correction applied as offset is supposed to be on top of drift region; however not at anode wire, so some inclination of drifting clusters due to Lorentz angle exists) are only !!
2073 // !!valid (in approximation) if tFS is close to the beginning of the drift region. !!
2074 // !!The slope however can be converted to a deflection length between electrode and cathode wire plane as it is clear that the drift region is sampled 20 times !!
2075 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2076
1d93b218 2077 // which formats should be chosen?
2078 AliTRDtrapAlu denomAlu;
2079 denomAlu.Init(20,8);
2080 AliTRDtrapAlu numAlu;
2081 numAlu.Init(20,8);
2082 // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
2083 // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
2084
2085 for (Int_t i = 0; i<4; i++) {
2086 if (mADC[i] == -1) continue;
2087
2088 Int_t num0 = (mN[i]*mXX[i]-mX[i]*mX[i]);
2089 if (num0 < 0) {
2090 denomAlu.AssignInt(-num0); // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2091 denomAlu.SetSign(-1);
2092 }
2093 else {
2094 denomAlu.AssignInt(num0);
2095 denomAlu.SetSign(1);
2096 }
2097
2098 Int_t num1 = mN[i]*mXY[i] - mX[i]*mY[i];
2099 if (num1 < 0) {
2100 numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2101 numAlu.SetSign(-1);
2102 }
2103 else {
2104 numAlu.AssignFormatted(num1);
2105 numAlu.SetSign(1);
2106 }
2107 numAlu = numAlu/denomAlu;
2108 mSlope[i] = numAlu.GetSignedValue();
2109
2110 Int_t num2 = mXX[i]*mY[i] - mX[i]*mXY[i];
2111
2112 if (num2 < 0) {
2113 numAlu.AssignFormatted(-num2);
2114 numAlu.SetSign(-1);
2115 }
2116 else {
2117 numAlu.AssignFormatted(num2);
2118 numAlu.SetSign(1);
2119 }
2120
2121 numAlu = numAlu/denomAlu;
2122
2123
2124 mOffset[i] = numAlu.GetSignedValue();
2125 numAlu.SetSign(1);
2126 denomAlu.SetSign(1);
2127
2128
2129 //numAlu.AssignInt(mADC[i]+1); // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)
2130 numAlu.AssignDouble((Double_t)mADC[i] + 1.5); // numAlu has enough pre-comma place for that; correct trafo, best values
2131 mOffset[i] = mOffset[i] + numAlu.GetValue(); // transform offset to a coord.system relative to chip; +1 to avoid neg. values
2132
96e6312d 2133 // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
1d93b218 2134 // reverse adc-counting to be again in line with the online mapping
2135 mADC[i] = fNADC - 4 - mADC[i]; // fNADC-4-mADC[i]: 0..17; remapping necessary;
2136 mADC[i] = mADC[i] + 2;
2137 // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2138 }
2139
2140 // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected;
2141 // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2142
2143 // transform parameters to the local coordinate-system of a stack (used by GTU)
2144 AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2145
2146 Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2147 //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2148
2149 // difference between width of inner and outer pads of a row is not accounted for;
2150
96e6312d 2151 Double_t magField = 0.4; // z-component of magnetic field in Tesla; adjust to current simulation!!; magnetic field can hardly be evaluated for the position of each mcm
2152 Double_t eCharge = 0.3; // unit charge in (GeV/c)/m*T
1d93b218 2153 Double_t ptMin = 2.3; // minimum transverse momentum (GeV/c); to be adjusted(?)
2154
2155 Double_t granularityOffset = 0.160; // granularity for offset in mm
2156 Double_t granularitySlope = 0.140; // granularity for slope in mm
2157
2158 // get the coordinates in SM-system; parameters:
2159
2160 Double_t zPos = (padPlane->GetRowPos(fRow))*10.0; // z-position of the MCM; fRow is counted on a chamber; SM consists of 5
2161 // zPos is position of pad-borders;
2162 Double_t zOffset = 0.0;
2163 if ( fRow == 0 || fRow == 15 ) {
2164 zOffset = padPlane->GetLengthOPad();
2165 }
2166 else {
2167 zOffset = padPlane->GetLengthIPad();
2168 }
2169 zOffset = (-1.0) * zOffset/2.0;
2170 // turn zPos to be z-coordinate at middle of pad-row
2171 zPos = zPos + zOffset;
2172
2173
2174 Double_t xPos = 0.0; // x-position of the upper border of the drift-chamber of actual layer
23200400 2175 Int_t icol = 0; // column-number of adc-channel
1d93b218 2176 Double_t yPos[4]; // y-position of the pad to which ADC is connected
23200400 2177 Double_t dx = 30.0; // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2178 Double_t freqSample = fFeeParam->GetSamplingFrequency(); // retrieve the sampling frequency (10.019750 MHz)
2179 Double_t vdrift = fCal->GetVdriftAverage(fChaId); // averaged drift velocity for this detector (1.500000 cm/us)
2180 Int_t nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
96e6312d 2181 Int_t nrOfAmplTimeBins = 2; // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
2182 Int_t nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
2183 if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed
1d93b218 2184 Double_t lorTan = fCal->GetOmegaTau(vdrift,magField); // tan of the Lorentz-angle for this detector; could be evaluated and set as a parameter for each mcm
23200400 2185 //Double_t lorAngle = 7.0; // Lorentz-angle in degrees
1d93b218 2186 Double_t tiltAngle = padPlane->GetTiltingAngle(); // sign-respecting tilting angle of pads in actual layer
2187 Double_t tiltTan = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2188 //Double_t lorTan = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2189
2190 Double_t alphaMax[4]; // maximum deflection from the direction to the primary vertex; granularity of hit pads
2191 Double_t slopeMin[4]; // local limits for the deflection
2192 Double_t slopeMax[4];
2193 Int_t mslopeMin[4]; // in granularity units; to be compared to mSlope[i]
2194 Int_t mslopeMax[4];
2195
2196
96e6312d 2197 // x coord. of upper side of drift chambers in local SM-system (in mm)
2198 // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2199 // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
1d93b218 2200 switch(fLayer)
2201 {
2202 case 0:
2203 xPos = 3003.0;
2204 break;
2205 case 1:
2206 xPos = 3129.0;
2207 break;
2208 case 2:
2209 xPos = 3255.0;
2210 break;
2211 case 3:
2212 xPos = 3381.0;
2213 break;
2214 case 4:
2215 xPos = 3507.0;
2216 break;
2217 case 5:
2218 xPos = 3633.0;
2219 break;
2220 }
2221
2222 // calculation of offset-correction n:
2223
2224 Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));
2225
2226 nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2227 if (nCorrectOffset < 0) {
2228 numAlu.AssignInt(-nCorrectOffset);
2229 numAlu.SetSign(-1);
2230 }
2231 else {
2232 numAlu.AssignInt(nCorrectOffset);
2233 numAlu.SetSign(1);
2234 }
96e6312d 2235 nCorrectOffset = numAlu.GetSignedValue();
2236
2237 // the Lorentz correction to the offset
2238 Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2239
2240
2241 lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2242
2243 if(lorCorrectOffset < 0) {
2244 numAlu.AssignDouble(-lorCorrectOffset);
2245 numAlu.SetSign(-1);
2246 }
2247 else{
2248 numAlu.AssignDouble(lorCorrectOffset);
2249 numAlu.SetSign(1);
2250 }
2251
2252 Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2253
2254
1d93b218 2255 Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2256
2257 // calculation of slope-correction
2258
2259 // this is only true for tracks coming (approx.) from primary vertex
23200400 2260 // everything is evaluated for a tracklet covering the whole drift chamber
1d93b218 2261 Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2262 // Double_t cCorrectSlope = zPos/xPos*dx*tiltTan/granularitySlope;
2263 // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2264
2265 if (cCorrectSlope < 0) {
2266 numAlu.AssignDouble(-cCorrectSlope);
2267 numAlu.SetSign(-1);
2268 }
2269 else {
2270 numAlu.AssignDouble(cCorrectSlope);
2271 numAlu.SetSign(1);
2272 }
2273 cCorrectSlope = numAlu.GetSignedValue();
2274
2275 // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2276 // different pad-width of outer pads of a pad-plane not taken into account
23200400 2277 // note that the fit was only done in the range tFS to tFE, however this range does not need to cover the whole drift region (neither start nor end of it)
2278 // however the tracklets are supposed to be a fit in the drift region thus the linear function is stretched to fit the drift region of 30 mm
2279
1d93b218 2280
23200400 2281 Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope; // >= 0.0
1d93b218 2282
2283 AliTRDtrapAlu correctAlu;
2284 correctAlu.Init(20,8);
2285
2286 AliTRDtrapAlu offsetAlu;
2287 offsetAlu.Init(13,0,-0x1000,0x0FFF); // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2288
2289 AliTRDtrapAlu slopeAlu;
2290 slopeAlu.Init(7,0,-0x40,0x3F); // 7 bit-word; 2-complement (1 sign-bit);
2291
2292 for (Int_t i = 0; i<4; i++) {
2293
2294 if (mADC[i] == -1) continue;
2295
2296 icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2297 yPos[i] = (padPlane->GetColPos(icol))*10.0;
2298
2299
2300 // offset:
2301
2302 correctAlu.AssignDouble(mCorrectOffset); // done because max. accuracy is 8 bit
2303 mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
96e6312d 2304 mOffset[i] = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset));
2305 //mOffset[i] = mOffset[i]*(-1); // adjust to direction of y-axes in online simulation
1d93b218 2306
2307 if (mOffset[i] < 0) {
2308 numAlu.AssignFormatted(-mOffset[i]);
2309 numAlu.SetSign(-1);
2310 }
2311 else {
2312 numAlu.AssignFormatted(mOffset[i]);
2313 numAlu.SetSign(1);
2314 }
2315
2316 offsetAlu = numAlu;
2317 mOffset[i] = offsetAlu.GetSignedValue();
2318
2319
2320 // slope:
2321
2322 correctAlu.AssignDouble(mCorrectSlope);
2323 mCorrectSlope = correctAlu.GetValueWhole();
2324
2325 mSlope[i] = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2326
2327 if (mSlope[i] < 0) {
2328 numAlu.AssignFormatted(-mSlope[i]);
2329 numAlu.SetSign(-1);
2330 }
2331 else {
2332 numAlu.AssignFormatted(mSlope[i]);
2333 numAlu.SetSign(1);
2334 }
2335
2336 slopeAlu = numAlu; // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2337 mSlope[i] = slopeAlu.GetSignedValue();
2338
2339 // local (LTU) limits for the deflection
2340 // ATan returns angles in radian
2341 alphaMax[i] = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2342 slopeMin[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2343 slopeMax[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2344
2345 if (slopeMin[i] < 0) {
2346 slopeAlu.AssignDouble(-slopeMin[i]);
2347 slopeAlu.SetSign(-1);
2348 }
2349 else {
2350 slopeAlu.AssignDouble(slopeMin[i]);
2351 slopeAlu.SetSign(1);
2352 }
2353 mslopeMin[i] = slopeAlu.GetSignedValue(); // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2354
2355 if (slopeMax[i] < 0) {
2356 slopeAlu.AssignDouble(-slopeMax[i]);
2357 slopeAlu.SetSign(-1);
2358 }
2359 else {
2360 slopeAlu.AssignDouble(slopeMax[i]);
2361 slopeAlu.SetSign(1);
2362 }
2363 mslopeMax[i] = slopeAlu.GetSignedValue();
2364 }
2365
2366 // suppress submission of tracks with low stiffness
2367 // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber)
2368
2369 // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2370 // up to now they are sorted according to maximum hit sum
2371 // is the sorting really done in the TRAP-chip?
2372
2373 Int_t order[4] = {-1,-1,-1,-1};
2374 Int_t wordnr = 0; // number of tracklet-words
2375
2376 for(Int_t j = 0; j < fMaxTracklets; j++) {
96e6312d 2377 //if( mADC[j] == -1) continue;
2378 if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
1d93b218 2379 wordnr++;
2380 if( wordnr-1 == 0) {
2381 order[0] = j;
2382 continue;
2383 }
2384 // wordnr-1>0, wordnr-1<4
2385 order[wordnr-1] = j;
2386 for( Int_t k = 0; k < wordnr-1; k++) {
2387 if( mOffset[j] < mOffset[order[k]] ) {
2388 for( Int_t l = wordnr-1; l > k; l-- ) {
2389 order[l] = order[l-1];
2390 }
2391 order[k] = j;
2392 break;
2393 }
2394
2395 }
2396 }
2397
2398 // fill the bit-words in ascending order and without gaps
2399 UInt_t bitWord[4] = {0,0,0,0}; // attention: unsigned int to have real 32 bits (no 2-complement)
2400 for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2401 //Bool_t rem1 = kTRUE;
2402
2403 Int_t i = order[j];
23200400 2404 //bit-word is 2-complement and therefore without sign
1d93b218 2405 bitWord[j] = 1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2406 //printf("\n");
23200400 2407 UInt_t shift = 0;
2408 UInt_t shift2 = 0;
2409
2410
1d93b218 2411
2412
2413 /*printf("mean charge: %d\n",mMeanCharge[i]);
2414 printf("row: %d\n",fRow);
2415 printf("slope: %d\n",mSlope[i]);
2416 printf("pad position: %d\n",mOffset[i]);
2417 printf("channel: %d\n",mADC[i]);*/
2418
2419 // electron probability (currently not implemented; the mean charge is just scaled)
23200400 2420 shift = (UInt_t)mMeanCharge[i];
1d93b218 2421 for(Int_t iBit = 0; iBit < 8; iBit++) {
2422 bitWord[j] = bitWord[j]<<1;
23200400 2423 bitWord[j] |= (shift>>(7-iBit))&1;
1d93b218 2424 //printf("0");
2425 }
23200400 2426
1d93b218 2427 // pad row
23200400 2428 shift = (UInt_t)fRow;
1d93b218 2429 for(Int_t iBit = 0; iBit < 4; iBit++) {
2430 bitWord[j] = bitWord[j]<<1;
23200400 2431 bitWord[j] |= (shift>>(3-iBit))&1;
1d93b218 2432 //printf("%d", (fRow>>(3-iBit))&1);
2433 }
2434
2435 // deflection length
1d93b218 2436 if(mSlope[i] < 0) {
23200400 2437 shift = (UInt_t)(-mSlope[i]);
2438 // shift2 is 2-complement of shift
2439 shift2 = 1;
2440 for(Int_t iBit = 1; iBit < 7; iBit++) {
2441 shift2 = shift2<<1;
5a67d089 2442 shift2 |= (1- (((shift)>>(6-iBit))&1) );
23200400 2443 //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2444 }
2445 shift2 = shift2 + 1;
2446 //printf("1");
2447 for(Int_t iBit = 0; iBit < 7; iBit++) {
2448 bitWord[j] = bitWord[j]<<1;
2449 bitWord[j] |= (shift2>>(6-iBit))&1;
2450 //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2451 }
1d93b218 2452 }
2453 else {
23200400 2454 shift = (UInt_t)(mSlope[i]);
1d93b218 2455 bitWord[j] = bitWord[j]<<1;
23200400 2456 bitWord[j] |= 0;
2457 //printf("0");
2458 for(Int_t iBit = 1; iBit < 7; iBit++) {
2459 bitWord[j] = bitWord[j]<<1;
2460 bitWord[j] |= (shift>>(6-iBit))&1;
2461 //printf("%d",(mSlope[i]>>(6-iBit))&1);
2462 }
1d93b218 2463 }
2464
2465 // pad position
1d93b218 2466 if(mOffset[i] < 0) {
23200400 2467 shift = (UInt_t)(-mOffset[i]);
2468 shift2 = 1;
1d93b218 2469 for(Int_t iBit = 1; iBit < 13; iBit++) {
23200400 2470 shift2 = shift2<<1;
5a67d089 2471 shift2 |= (1-(((shift)>>(12-iBit))&1));
1d93b218 2472 //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2473 }
23200400 2474 shift2 = shift2 + 1;
2475 //printf("1");
2476 for(Int_t iBit = 0; iBit < 13; iBit++) {
2477 bitWord[j] = bitWord[j]<<1;
2478 bitWord[j] |= (shift2>>(12-iBit))&1;
2479 //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2480 }
1d93b218 2481 }
2482 else {
23200400 2483 shift = (UInt_t)mOffset[i];
2484 bitWord[j] = bitWord[j]<<1;
2485 bitWord[j] |= 0;
2486 //printf("0");
2487 for(Int_t iBit = 1; iBit < 13; iBit++) {
2488 bitWord[j] = bitWord[j]<<1;
2489 bitWord[j] |= (shift>>(12-iBit))&1;
2490 //printf("%d",(mOffset[i]>>(12-iBit))&1);
2491 }
1d93b218 2492 }
2493
2494
23200400 2495
1d93b218 2496 //printf("bitWord: %u\n",bitWord[j]);
2497 //printf("adc: %d\n",mADC[i]);
2498 fMCMT[j] = bitWord[j];
2499 }
2500
2501 //printf("\n");
2502
2503
2504 delete [] qsum;
2505 delete [] ieffped;
2506
2507 delete [] X;
2508 delete [] XX;
2509 delete [] Y;
2510 delete [] YY;
2511 delete [] XY;
2512 delete [] N;
2513 delete [] QT0;
2514 delete [] QT1;
2515
2516 delete [] hitSum;
2517 delete [] selectPair;
2518
2519 delete padPlane;
2520
96e6312d 2521//if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2522
4cc89512 2523 if (!fFeeParam->GetMCTrackletOutput())
2524 return;
1d93b218 2525
1d93b218 2526 AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2527 AliLog::SetFileOutput("../log/tracklet.log");
2528
1d93b218 2529 // testing for wordnr in order to speed up the simulation
52c19022 2530 if (wordnr == 0)
1d93b218 2531 return;
1d93b218 2532
4cc89512 2533 UInt_t *trackletWord = new UInt_t[fMaxTracklets];
2534 Int_t *adcChannel = new Int_t[fMaxTracklets];
2535 Int_t *trackRef = new Int_t[fMaxTracklets];
2536
2537 Int_t u = 0;
2538
2539 AliTRDdigitsManager *digman = new AliTRDdigitsManager();
2540 digman->ReadDigits(gAlice->GetRunLoader()->GetLoader("TRDLoader")->TreeD());
2541 digman->SetUseDictionaries(kTRUE);
2542 AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
1d93b218 2543
2544 for (Int_t j = 0; j < fMaxTracklets; j++) {
2545 Int_t i = order[j];
2546 trackletWord[j] = 0;
2547 adcChannel[j] = -1;
2548 if (bitWord[j]!=0) {
2549 trackletWord[u] = bitWord[j];
2550 adcChannel[u] = mADC[i]; // mapping onto the original adc-array to be in line with the digits-adc-ordering (21 channels in total on 1 mcm, 18 belonging to pads); mADC[i] should be >-1 in case bitWord[i]>0
4cc89512 2551
2552// Finding label of MC track
2553 TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
2554 Int_t track[3];
2555 Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
2556 Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
2557 Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
2558 Int_t det = 30 * fSector + 6 * fStack + fLayer;
2559 for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
2560 track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2561 track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2562 track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2563 hTrkRef->Fill(track[0]);
2564 if (track[1] != track[0] && track[1] != -1)
2565 hTrkRef->Fill(track[1]);
2566 if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2567 hTrkRef->Fill(track[2]);
2568 if (padcol_ngb >= 0) {
2569 track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2570 track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2571 track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2572 hTrkRef->Fill(track[0]);
2573 if (track[1] != track[0] && track[1] != -1)
2574 hTrkRef->Fill(track[1]);
2575 if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2576 hTrkRef->Fill(track[2]);
2577 }
2578 }
2579 trackRef[u] = hTrkRef->GetMaximumBin() - 1;
2580 delete hTrkRef;
1d93b218 2581 u = u + 1;
2582 }
2583 }
1d93b218 2584
52c19022 2585 AliDataLoader *dl = gAlice->GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
2586 if (!dl) {
2587 AliError("Could not get the tracklets data loader!");
1d93b218 2588 }
52c19022 2589 else {
2590 TTree *trackletTree = dl->Tree();
2591 if (!trackletTree)
2592 dl->MakeTree();
2593 trackletTree = dl->Tree();
1d93b218 2594
52c19022 2595 AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM();
2596 TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
2597 if (!trkbranch)
2598 trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
2599 trkbranch->SetAddress(&trkl);
1d93b218 2600
52c19022 2601 for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
2602 if (trackletWord[iTracklet] == 0)
2603 continue;
2604 trkl->SetTrackletWord(trackletWord[iTracklet]);
2605 trkl->SetDetector(30*fSector + 6*fStack + fLayer);
2606 trkl->SetROB(fRobPos);
2607 trkl->SetMCM(fMcmPos);
4cc89512 2608 trkl->SetLabel(trackRef[iTracklet]);
52c19022 2609 trackletTree->Fill();
52c19022 2610 }
2611 delete trkl;
2612 dl->WriteData("OVERWRITE");
2613 }
23200400 2614
4cc89512 2615 delete [] trackletWord;
2616 delete [] adcChannel;
2617 delete [] trackRef;
2618 delete digman;
1d93b218 2619
2620 // to be done:
2621 // error measure for quality of fit (not necessarily needed for the trigger)
2622 // cluster quality threshold (not yet set)
2623 // electron probability
1d93b218 2624}
2625