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