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