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