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