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