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