Making CINT happy
[u/mrichter/AliRoot.git] / TRD / AliTRDdigitizer.cxx
CommitLineData
f7336fa3 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
16/*
17$Log$
3583ac1a 18Revision 1.40 2002/10/14 14:57:43 hristov
19Merging the VirtualMC branch to the main development branch (HEAD)
20
b9d0a01d 21Revision 1.33.6.3 2002/10/11 07:26:37 hristov
22Updating VirtualMC to v3-09-02
23
24Revision 1.39 2002/10/08 20:46:12 cblume
25Do coupling factors before noise is applied
26
0c24ba98 27Revision 1.38 2002/04/30 08:30:40 cblume
28gAlice now only read by AliRunDigitizer. Therefore it is just deleted in AliTRDmerge.C
29
c57e2264 30Revision 1.37 2002/04/29 11:50:47 cblume
31Change initialization of gAlice in the merging case
32
cec4059b 33Revision 1.36 2002/04/12 12:13:23 cblume
34Add Jiris changes
35
8e64dd77 36Revision 1.35 2002/03/28 14:59:07 cblume
37Coding conventions
38
0a29d0f1 39Revision 1.34 2002/03/25 20:00:44 cblume
40Introduce parameter class and regions of interest for merging
41
4487dad0 42Revision 1.33 2002/02/12 17:32:03 cblume
43Rearrange the deleting of the list of sdigitsmanager
44
b1113c6b 45Revision 1.32 2002/02/12 16:07:21 cblume
46Add new constructor
47
bfc40adc 48Revision 1.31 2002/02/11 14:27:11 cblume
49New pad plane design, new TRF+PRF, tail cancelation, cross talk
50
e23fbb27 51Revision 1.30 2001/11/19 08:44:08 cblume
52Fix bugs reported by Rene
53
a1c3aded 54Revision 1.29 2001/11/14 19:44:25 hristov
55Numeric const casted (Alpha)
56
20c69486 57Revision 1.28 2001/11/14 16:35:58 cblume
58Inherits now from AliDetector
59
85cbec76 60Revision 1.27 2001/11/14 10:50:45 cblume
61Changes in digits IO. Add merging of summable digits
62
abaf1f1d 63Revision 1.26 2001/11/06 17:19:41 cblume
64Add detailed geometry and simple simulator
65
16bf9884 66Revision 1.25 2001/06/27 09:54:44 cblume
67Moved fField initialization to InitDetector()
68
88fd7ccb 69Revision 1.24 2001/05/21 16:45:47 hristov
70Last minute changes (C.Blume)
71
db30bf0f 72Revision 1.23 2001/05/07 08:04:48 cblume
73New TRF and PRF. Speedup of the code. Digits from amplification region included
74
872a7aba 75Revision 1.22 2001/03/30 14:40:14 cblume
76Update of the digitization parameter
77
a3c76cdc 78Revision 1.21 2001/03/13 09:30:35 cblume
79Update of digitization. Moved digit branch definition to AliTRD
80
6244debe 81Revision 1.20 2001/02/25 20:19:00 hristov
82Minor correction: loop variable declared only once for HP, Sun
83
c3a4830f 84Revision 1.19 2001/02/14 18:22:26 cblume
85Change in the geometry of the padplane
86
71d9fa7b 87Revision 1.18 2001/01/26 19:56:57 hristov
88Major upgrade of AliRoot code
89
2ab0c725 90Revision 1.17 2000/12/08 12:53:27 cblume
91Change in Copy() function for HP-compiler
92
1948ba0c 93Revision 1.16 2000/12/07 12:20:46 cblume
94Go back to array compression. Use sampled PRF to speed up digitization
95
e153aaf6 96Revision 1.15 2000/11/23 14:34:08 cblume
97Fixed bug in expansion routine of arrays (initialize buffers properly)
98
259b9e4b 99Revision 1.14 2000/11/20 08:54:44 cblume
100Switch off compression as default
101
c1e4b257 102Revision 1.13 2000/11/10 14:57:52 cblume
103Changes in the geometry constants for the DEC compiler
104
dd56b762 105Revision 1.12 2000/11/01 14:53:20 cblume
106Merge with TRD-develop
107
793ff80c 108Revision 1.1.4.9 2000/10/26 17:00:22 cblume
109Fixed bug in CheckDetector()
110
111Revision 1.1.4.8 2000/10/23 13:41:35 cblume
112Added protection against Log(0) in the gas gain calulation
113
114Revision 1.1.4.7 2000/10/17 02:27:34 cblume
115Get rid of global constants
116
117Revision 1.1.4.6 2000/10/16 01:16:53 cblume
118Changed timebin 0 to be the one closest to the readout
119
120Revision 1.1.4.5 2000/10/15 23:34:29 cblume
121Faster version of the digitizer
122
123Revision 1.1.4.4 2000/10/06 16:49:46 cblume
124Made Getters const
125
126Revision 1.1.4.3 2000/10/04 16:34:58 cblume
127Replace include files by forward declarations
128
129Revision 1.1.4.2 2000/09/22 14:41:10 cblume
130Bug fix in PRF. Included time response. New structure
131
eda4336d 132Revision 1.10 2000/10/05 07:27:53 cblume
133Changes in the header-files by FCA
134
6798b56e 135Revision 1.9 2000/10/02 21:28:19 fca
136Removal of useless dependecies via forward declarations
137
94de3818 138Revision 1.8 2000/06/09 11:10:07 cblume
139Compiler warnings and coding conventions, next round
140
dd9a6ee3 141Revision 1.7 2000/06/08 18:32:58 cblume
142Make code compliant to coding conventions
143
8230f242 144Revision 1.6 2000/06/07 16:27:32 cblume
145Try to remove compiler warnings on Sun and HP
146
9d0b222b 147Revision 1.5 2000/05/09 16:38:57 cblume
148Removed PadResponse(). Merge problem
149
f0a7bf65 150Revision 1.4 2000/05/08 15:53:45 cblume
151Resolved merge conflict
152
da581aea 153Revision 1.3 2000/04/28 14:49:27 cblume
154Only one declaration of iDict in MakeDigits()
155
156Revision 1.1.4.1 2000/05/08 14:42:04 cblume
157Introduced AliTRDdigitsManager
28329a48 158
1befd3b2 159Revision 1.1 2000/02/28 19:00:13 cblume
160Add new TRD classes
161
f7336fa3 162*/
163
164///////////////////////////////////////////////////////////////////////////////
165// //
166// Creates and handles digits from TRD hits //
abaf1f1d 167// Author: C. Blume (C.Blume@gsi.de) //
f7336fa3 168// //
169// The following effects are included: //
170// - Diffusion //
171// - ExB effects //
172// - Gas gain including fluctuations //
173// - Pad-response (simple Gaussian approximation) //
abaf1f1d 174// - Time-response //
f7336fa3 175// - Electronics noise //
176// - Electronics gain //
177// - Digitization //
178// - ADC threshold //
179// The corresponding parameter can be adjusted via the various //
180// Set-functions. If these parameters are not explicitly set, default //
181// values are used (see Init-function). //
abaf1f1d 182// As an example on how to use this class to produce digits from hits //
183// have a look at the macro hits2digits.C //
184// The production of summable digits is demonstrated in hits2sdigits.C //
185// and the subsequent conversion of the s-digits into normal digits is //
186// explained in sdigits2digits.C. //
f7336fa3 187// //
188///////////////////////////////////////////////////////////////////////////////
189
6798b56e 190#include <stdlib.h>
191
f7336fa3 192#include <TMath.h>
193#include <TVector.h>
194#include <TRandom.h>
94de3818 195#include <TROOT.h>
196#include <TTree.h>
793ff80c 197#include <TFile.h>
198#include <TF1.h>
abaf1f1d 199#include <TList.h>
85cbec76 200#include <TTask.h>
793ff80c 201
202#include "AliRun.h"
db30bf0f 203#include "AliMagF.h"
85cbec76 204#include "AliRunDigitizer.h"
f7336fa3 205
206#include "AliTRD.h"
793ff80c 207#include "AliTRDhit.h"
f7336fa3 208#include "AliTRDdigitizer.h"
da581aea 209#include "AliTRDdataArrayI.h"
210#include "AliTRDdataArrayF.h"
793ff80c 211#include "AliTRDsegmentArray.h"
da581aea 212#include "AliTRDdigitsManager.h"
793ff80c 213#include "AliTRDgeometry.h"
4487dad0 214#include "AliTRDparameter.h"
f7336fa3 215
216ClassImp(AliTRDdigitizer)
217
218//_____________________________________________________________________________
85cbec76 219AliTRDdigitizer::AliTRDdigitizer()
f7336fa3 220{
221 //
222 // AliTRDdigitizer default constructor
223 //
224
4487dad0 225 fInputFile = 0;
226 fDigitsManager = 0;
227 fSDigitsManagerList = 0;
228 fSDigitsManager = 0;
229 fTRD = 0;
230 fGeo = 0;
231 fPar = 0;
e23fbb27 232 fMasks = 0;
abaf1f1d 233 fEvent = 0;
abaf1f1d 234 fCompress = kTRUE;
e23fbb27 235 fDebug = 0;
abaf1f1d 236 fSDigits = kFALSE;
237 fSDigitsScale = 0.0;
4487dad0 238 fMergeSignalOnly = kFALSE;
0c24ba98 239 fSimpleSim = kFALSE;
240 fSimpleDet = 0;
241
f7336fa3 242}
243
244//_____________________________________________________________________________
245AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
85cbec76 246 :AliDigitizer(name,title)
f7336fa3 247{
248 //
85cbec76 249 // AliTRDdigitizer constructor
250 //
251
4487dad0 252 fInputFile = 0;
253 fDigitsManager = 0;
254 fSDigitsManagerList = 0;
255 fSDigitsManager = 0;
256 fTRD = 0;
257 fGeo = 0;
258 fPar = 0;
e23fbb27 259 fMasks = 0;
85cbec76 260 fEvent = 0;
85cbec76 261 fCompress = kTRUE;
e23fbb27 262 fDebug = 0;
85cbec76 263 fSDigits = kFALSE;
4487dad0 264 fMergeSignalOnly = kFALSE;
0c24ba98 265 fSimpleSim = kFALSE;
266 fSimpleDet = 0;
267
4487dad0 268 // For the summable digits
269 fSDigitsScale = 100.;
85cbec76 270
271}
272
273//_____________________________________________________________________________
274AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
275 , const Text_t *name, const Text_t *title)
276 :AliDigitizer(manager,name,title)
277{
278 //
279 // AliTRDdigitizer constructor
f7336fa3 280 //
281
4487dad0 282 fInputFile = 0;
283 fDigitsManager = 0;
284 fSDigitsManagerList = 0;
285 fSDigitsManager = 0;
286 fTRD = 0;
287 fGeo = 0;
288 fPar = 0;
bfc40adc 289 fMasks = 0;
bfc40adc 290 fEvent = 0;
bfc40adc 291 fCompress = kTRUE;
292 fDebug = 0;
293 fSDigits = kFALSE;
4487dad0 294 fMergeSignalOnly = kFALSE;
0c24ba98 295 fSimpleSim = kFALSE;
296 fSimpleDet = 0;
297
4487dad0 298 // For the summable digits
299 fSDigitsScale = 100.;
bfc40adc 300
301}
302
303//_____________________________________________________________________________
304AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
305 :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
306{
307 //
308 // AliTRDdigitizer constructor
309 //
310
4487dad0 311 fInputFile = 0;
312 fDigitsManager = 0;
313 fSDigitsManagerList = 0;
314 fSDigitsManager = 0;
315 fTRD = 0;
316 fGeo = 0;
317 fPar = 0;
e23fbb27 318 fMasks = 0;
abaf1f1d 319 fEvent = 0;
abaf1f1d 320 fCompress = kTRUE;
e23fbb27 321 fDebug = 0;
abaf1f1d 322 fSDigits = kFALSE;
4487dad0 323 fMergeSignalOnly = kFALSE;
0c24ba98 324 fSimpleSim = kFALSE;
325 fSimpleDet = 0;
793ff80c 326
4487dad0 327 // For the summable digits
328 fSDigitsScale = 100.;
f7336fa3 329
330}
331
332//_____________________________________________________________________________
dd9a6ee3 333AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
8230f242 334{
335 //
336 // AliTRDdigitizer copy constructor
337 //
338
dd9a6ee3 339 ((AliTRDdigitizer &) d).Copy(*this);
8230f242 340
341}
342
343//_____________________________________________________________________________
f7336fa3 344AliTRDdigitizer::~AliTRDdigitizer()
345{
8230f242 346 //
347 // AliTRDdigitizer destructor
348 //
f7336fa3 349
350 if (fInputFile) {
351 fInputFile->Close();
352 delete fInputFile;
4487dad0 353 fInputFile = 0;
abaf1f1d 354 }
355
356 if (fDigitsManager) {
357 delete fDigitsManager;
4487dad0 358 fDigitsManager = 0;
f7336fa3 359 }
360
abaf1f1d 361 if (fSDigitsManager) {
362 delete fSDigitsManager;
4487dad0 363 fSDigitsManager = 0;
f7336fa3 364 }
365
abaf1f1d 366 if (fSDigitsManagerList) {
abaf1f1d 367 delete fSDigitsManagerList;
4487dad0 368 fSDigitsManagerList = 0;
e23fbb27 369 }
370
371 if (fMasks) {
372 delete [] fMasks;
373 fMasks = 0;
374 }
375
f7336fa3 376}
377
378//_____________________________________________________________________________
dd9a6ee3 379AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
380{
381 //
382 // Assignment operator
383 //
384
385 if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
386 return *this;
387
388}
389
390//_____________________________________________________________________________
391void AliTRDdigitizer::Copy(TObject &d)
8230f242 392{
393 //
394 // Copy function
395 //
396
4487dad0 397 ((AliTRDdigitizer &) d).fInputFile = 0;
398 ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
399 ((AliTRDdigitizer &) d).fSDigitsManager = 0;
400 ((AliTRDdigitizer &) d).fDigitsManager = 0;
401 ((AliTRDdigitizer &) d).fTRD = 0;
402 ((AliTRDdigitizer &) d).fGeo = 0;
e23fbb27 403 ((AliTRDdigitizer &) d).fMasks = 0;
abaf1f1d 404 ((AliTRDdigitizer &) d).fEvent = 0;
4487dad0 405 ((AliTRDdigitizer &) d).fPar = 0;
abaf1f1d 406 ((AliTRDdigitizer &) d).fCompress = fCompress;
e23fbb27 407 ((AliTRDdigitizer &) d).fDebug = fDebug ;
abaf1f1d 408 ((AliTRDdigitizer &) d).fSDigits = fSDigits;
409 ((AliTRDdigitizer &) d).fSDigitsScale = fSDigitsScale;
4487dad0 410 ((AliTRDdigitizer &) d).fMergeSignalOnly = fMergeSignalOnly;
0c24ba98 411 ((AliTRDdigitizer &) d).fSimpleSim = fSimpleSim;
412 ((AliTRDdigitizer &) d).fSimpleDet = fSimpleDet;
413
e23fbb27 414}
415
416//_____________________________________________________________________________
e23fbb27 417void AliTRDdigitizer::Exec(Option_t* option)
418{
419 //
420 // Executes the merging
421 //
422
423 Int_t iInput;
424
425 AliTRDdigitsManager *sdigitsManager;
426
427 TString optionString = option;
428 if (optionString.Contains("deb")) {
429 fDebug = 1;
430 if (optionString.Contains("2")) {
431 fDebug = 2;
432 }
433 printf("<AliTRDdigitizer::Exec> ");
434 printf("Called with debug option %d\n",fDebug);
435 }
436
cec4059b 437 // The AliRoot file is already connected by the manager
4487dad0 438 if (gAlice) {
439 if (fDebug > 0) {
440 printf("<AliTRDdigitizer::Exec> ");
441 printf("AliRun object found on file.\n");
442 }
443 }
444 else {
445 printf("<AliTRDdigitizer::Exec> ");
446 printf("Could not find AliRun object.\n");
cec4059b 447 return;
4487dad0 448 }
cec4059b 449
e23fbb27 450 Int_t nInput = fManager->GetNinputs();
451 fMasks = new Int_t[nInput];
452 for (iInput = 0; iInput < nInput; iInput++) {
453 fMasks[iInput] = fManager->GetMask(iInput);
454 }
e23fbb27 455
456 // Initialization
457 InitDetector();
458
459 for (iInput = 0; iInput < nInput; iInput++) {
460
461 if (fDebug > 0) {
462 printf("<AliTRDdigitizer::Exec> ");
463 printf("Add input stream %d\n",iInput);
464 }
465
c57e2264 466 // check if the input tree exists
467 if (!fManager->GetInputTreeTRDS(iInput)) {
468 printf("<AliTRDdigitizer::Exec> ");
469 printf("Input stream %d does not exist\n",iInput);
470 return;
471 }
472
e23fbb27 473 // Read the s-digits via digits manager
474 sdigitsManager = new AliTRDdigitsManager();
475 sdigitsManager->SetDebug(fDebug);
476 sdigitsManager->SetSDigits(kTRUE);
477 sdigitsManager->ReadDigits(fManager->GetInputTreeTRDS(iInput));
478
479 // Add the s-digits to the input list
480 AddSDigitsManager(sdigitsManager);
481
482 }
483
484 // Convert the s-digits to normal digits
485 if (fDebug > 0) {
486 printf("<AliTRDdigitizer::Exec> ");
487 printf("Do the conversion\n");
488 }
489 SDigits2Digits();
490
491 // Store the digits
492 if (fDebug > 0) {
493 printf("<AliTRDdigitizer::Exec> ");
494 printf("Write the digits\n");
495 }
496 fDigitsManager->MakeBranch(fManager->GetTreeDTRD());
497 fDigitsManager->WriteDigits();
498 if (fDebug > 0) {
499 printf("<AliTRDdigitizer::Exec> ");
500 printf("Done\n");
501 }
502
b1113c6b 503 DeleteSDigitsManager();
504
e23fbb27 505}
506
507//_____________________________________________________________________________
e23fbb27 508Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
f7336fa3 509{
510 //
511 // Opens a ROOT-file with TRD-hits and reads in the hit-tree
512 //
513
514 // Connect the AliRoot file containing Geometry, Kine, and Hits
4487dad0 515 fInputFile = (TFile *) gROOT->GetListOfFiles()->FindObject(file);
f7336fa3 516 if (!fInputFile) {
e23fbb27 517 if (fDebug > 0) {
518 printf("<AliTRDdigitizer::Open> ");
519 printf("Open the AliROOT-file %s.\n",file);
abaf1f1d 520 }
e23fbb27 521 fInputFile = new TFile(file,"UPDATE");
f7336fa3 522 }
523 else {
e23fbb27 524 if (fDebug > 0) {
525 printf("<AliTRDdigitizer::Open> ");
526 printf("%s is already open.\n",file);
abaf1f1d 527 }
f7336fa3 528 }
529
4487dad0 530 gAlice = (AliRun *) fInputFile->Get("gAlice");
da581aea 531 if (gAlice) {
e23fbb27 532 if (fDebug > 0) {
533 printf("<AliTRDdigitizer::Open> ");
abaf1f1d 534 printf("AliRun object found on file.\n");
535 }
da581aea 536 }
537 else {
e23fbb27 538 printf("<AliTRDdigitizer::Open> ");
da581aea 539 printf("Could not find AliRun object.\n");
540 return kFALSE;
541 }
f7336fa3 542
543 fEvent = nEvent;
544
545 // Import the Trees for the event nEvent in the file
546 Int_t nparticles = gAlice->GetEvent(fEvent);
547 if (nparticles <= 0) {
e23fbb27 548 printf("<AliTRDdigitizer::Open> ");
f7336fa3 549 printf("No entries in the trees for event %d.\n",fEvent);
550 return kFALSE;
551 }
552
abaf1f1d 553 if (InitDetector()) {
554 return MakeBranch();
555 }
556 else {
557 return kFALSE;
558 }
793ff80c 559
560}
561
562//_____________________________________________________________________________
563Bool_t AliTRDdigitizer::InitDetector()
564{
565 //
566 // Sets the pointer to the TRD detector and the geometry
567 //
568
dd9a6ee3 569 // Get the pointer to the detector class and check for version 1
4487dad0 570 fTRD = (AliTRD *) gAlice->GetDetector("TRD");
cec4059b 571 if (!fTRD) {
572 printf("<AliTRDdigitizer::InitDetector> ");
573 printf("No TRD module found\n");
574 exit(1);
575 }
dd9a6ee3 576 if (fTRD->IsVersion() != 1) {
e23fbb27 577 printf("<AliTRDdigitizer::InitDetector> ");
dd9a6ee3 578 printf("TRD must be version 1 (slow simulator).\n");
579 exit(1);
580 }
581
582 // Get the geometry
583 fGeo = fTRD->GetGeometry();
e23fbb27 584 if (fDebug > 0) {
585 printf("<AliTRDdigitizer::InitDetector> ");
abaf1f1d 586 printf("Geometry version %d\n",fGeo->IsVersion());
587 }
dd9a6ee3 588
abaf1f1d 589 // Create a digits manager
590 fDigitsManager = new AliTRDdigitsManager();
591 fDigitsManager->SetSDigits(fSDigits);
592 fDigitsManager->CreateArrays();
593 fDigitsManager->SetEvent(fEvent);
e23fbb27 594 fDigitsManager->SetDebug(fDebug);
abaf1f1d 595
596 // The list for the input s-digits manager to be merged
597 fSDigitsManagerList = new TList();
598
4487dad0 599 return kTRUE;
f7336fa3 600
601}
602
603//_____________________________________________________________________________
0a29d0f1 604Bool_t AliTRDdigitizer::MakeBranch(const Char_t *file) const
6244debe 605{
abaf1f1d 606 //
607 // Create the branches for the digits array
6244debe 608 //
609
abaf1f1d 610 return fDigitsManager->MakeBranch(file);
6244debe 611
612}
613
614//_____________________________________________________________________________
f7336fa3 615Bool_t AliTRDdigitizer::MakeDigits()
616{
617 //
872a7aba 618 // Creates digits.
f7336fa3 619 //
620
f7336fa3 621 ///////////////////////////////////////////////////////////////
622 // Parameter
623 ///////////////////////////////////////////////////////////////
624
625 // Converts number of electrons to fC
872a7aba 626 const Double_t kEl2fC = 1.602E-19 * 1.0E15;
f7336fa3 627
628 ///////////////////////////////////////////////////////////////
629
793ff80c 630 // Number of pads included in the pad response
631 const Int_t kNpad = 3;
632
633 // Number of track dictionary arrays
dd56b762 634 const Int_t kNDict = AliTRDdigitsManager::kNDict;
793ff80c 635
872a7aba 636 // Half the width of the amplification region
637 const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
638
c3a4830f 639 Int_t iRow, iCol, iTime, iPad;
71d9fa7b 640 Int_t iDict = 0;
793ff80c 641 Int_t nBytes = 0;
f7336fa3 642
643 Int_t totalSizeDigits = 0;
644 Int_t totalSizeDict0 = 0;
645 Int_t totalSizeDict1 = 0;
646 Int_t totalSizeDict2 = 0;
647
872a7aba 648 Int_t timeTRDbeg = 0;
649 Int_t timeTRDend = 1;
650
651 Float_t pos[3];
652 Float_t rot[3];
653 Float_t xyz[3];
654 Float_t padSignal[kNpad];
655 Float_t signalOld[kNpad];
656
793ff80c 657 AliTRDdataArrayF *signals = 0;
658 AliTRDdataArrayI *digits = 0;
8230f242 659 AliTRDdataArrayI *dictionary[kNDict];
da581aea 660
4487dad0 661 // Create a default parameter class if none is defined
662 if (!fPar) {
663 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
664 if (fDebug > 0) {
665 printf("<AliTRDdigitizer::MakeDigits> ");
666 printf("Create the default parameter object\n");
667 }
668 }
669
793ff80c 670 // Create a container for the amplitudes
671 AliTRDsegmentArray *signalsArray
e23fbb27 672 = new AliTRDsegmentArray("AliTRDdataArrayF"
673 ,AliTRDgeometry::Ndet());
793ff80c 674
4487dad0 675 if (fPar->TRFOn()) {
676 timeTRDbeg = ((Int_t) (-fPar->GetTRFlo() / fPar->GetTimeBinSize())) - 1;
677 timeTRDend = ((Int_t) ( fPar->GetTRFhi() / fPar->GetTimeBinSize())) - 1;
e23fbb27 678 if (fDebug > 0) {
679 printf("<AliTRDdigitizer::MakeDigits> ");
abaf1f1d 680 printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
681 }
da581aea 682 }
683
4487dad0 684 Float_t elAttachProp = fPar->GetElAttachProp() / 100.;
793ff80c 685
872a7aba 686 if (!fGeo) {
e23fbb27 687 printf("<AliTRDdigitizer::MakeDigits> ");
872a7aba 688 printf("No geometry defined\n");
689 return kFALSE;
690 }
691
e23fbb27 692 if (fDebug > 0) {
693 printf("<AliTRDdigitizer::MakeDigits> ");
abaf1f1d 694 printf("Start creating digits.\n");
695 }
872a7aba 696
793ff80c 697 // Get the pointer to the hit tree
0a29d0f1 698 TTree *hitTree = gAlice->TreeH();
793ff80c 699
700 // Get the number of entries in the hit tree
701 // (Number of primary particles creating a hit somewhere)
0c24ba98 702 Int_t nTrack = 1;
703 if (!fSimpleSim) {
704 nTrack = (Int_t) hitTree->GetEntries();
705 if (fDebug > 0) {
706 printf("<AliTRDdigitizer::MakeDigits> ");
707 printf("Found %d primary particles\n",nTrack);
708 }
709 }
793ff80c 710
711 Int_t detectorOld = -1;
712 Int_t countHits = 0;
713
714 // Loop through all entries in the tree
715 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
716
0c24ba98 717 if (!fSimpleSim) {
718 gAlice->ResetHits();
719 nBytes += hitTree->GetEvent(iTrack);
720 }
793ff80c 721
e23fbb27 722 // Loop through the TRD hits
723 Int_t iHit = 0;
724 AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
725 while (hit) {
726
793ff80c 727 countHits++;
e23fbb27 728 iHit++;
729
730 pos[0] = hit->X();
731 pos[1] = hit->Y();
732 pos[2] = hit->Z();
733 Float_t q = hit->GetCharge();
734 Int_t track = hit->Track();
735 Int_t detector = hit->GetDetector();
736 Int_t plane = fGeo->GetPlane(detector);
737 Int_t sector = fGeo->GetSector(detector);
738 Int_t chamber = fGeo->GetChamber(detector);
4487dad0 739 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
740 Int_t nColMax = fPar->GetColMax(plane);
741 Int_t nTimeMax = fPar->GetTimeMax();
742 Int_t nTimeBefore = fPar->GetTimeBefore();
743 Int_t nTimeAfter = fPar->GetTimeAfter();
744 Int_t nTimeTotal = fPar->GetTimeTotal();
745 Float_t row0 = fPar->GetRow0(plane,chamber,sector);
746 Float_t col0 = fPar->GetCol0(plane);
747 Float_t time0 = fPar->GetTime0(plane);
748 Float_t rowPadSize = fPar->GetRowPadSize(plane,chamber,sector);
749 Float_t colPadSize = fPar->GetColPadSize(plane);
750 Float_t timeBinSize = fPar->GetTimeBinSize();
872a7aba 751 Float_t divideRow = 1.0 / rowPadSize;
752 Float_t divideCol = 1.0 / colPadSize;
753 Float_t divideTime = 1.0 / timeBinSize;
793ff80c 754
e23fbb27 755 if (fDebug > 1) {
793ff80c 756 printf("Analyze hit no. %d ",iHit);
757 printf("-----------------------------------------------------------\n");
758 hit->Dump();
759 printf("plane = %d, sector = %d, chamber = %d\n"
760 ,plane,sector,chamber);
761 printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n"
762 ,nRowMax,nColMax,nTimeMax);
872a7aba 763 printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
764 ,nTimeBefore,nTimeAfter,nTimeTotal);
793ff80c 765 printf("row0 = %f, col0 = %f, time0 = %f\n"
766 ,row0,col0,time0);
c1e4b257 767 printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
768 ,rowPadSize,colPadSize,timeBinSize);
793ff80c 769 }
770
e23fbb27 771 // Don't analyze test hits and switched off detectors
772 if ((CheckDetector(plane,chamber,sector)) &&
773 (((Int_t) q) != 0)) {
dd56b762 774
e23fbb27 775 if (detector != detectorOld) {
e153aaf6 776
e23fbb27 777 if (fDebug > 1) {
778 printf("<AliTRDdigitizer::MakeDigits> ");
779 printf("Get new container. New det = %d, Old det = %d\n"
780 ,detector,detectorOld);
dd9a6ee3 781 }
e23fbb27 782 // Compress the old one if enabled
783 if ((fCompress) && (detectorOld > -1)) {
784 if (fDebug > 1) {
785 printf("<AliTRDdigitizer::MakeDigits> ");
786 printf("Compress the old container ...");
787 }
788 signals->Compress(1,0);
789 for (iDict = 0; iDict < kNDict; iDict++) {
790 dictionary[iDict]->Compress(1,0);
c1e4b257 791 }
e23fbb27 792 if (fDebug > 1) printf("done\n");
793ff80c 793 }
e23fbb27 794 // Get the new container
795 signals = (AliTRDdataArrayF *) signalsArray->At(detector);
796 if (signals->GetNtime() == 0) {
797 // Allocate a new one if not yet existing
798 if (fDebug > 1) {
799 printf("<AliTRDdigitizer::MakeDigits> ");
800 printf("Allocate a new container ... ");
801 }
802 signals->Allocate(nRowMax,nColMax,nTimeTotal);
793ff80c 803 }
0c24ba98 804 else if (fSimpleSim) {
805 // Clear an old one for the simple simulation
806 if (fDebug > 1) {
807 printf("<AliTRDdigitizer::MakeDigits> ");
808 printf("Clear a old container ... ");
809 }
810 signals->Clear();
811 }
793ff80c 812 else {
e23fbb27 813 // Expand an existing one
814 if (fCompress) {
815 if (fDebug > 1) {
816 printf("<AliTRDdigitizer::MakeDigits> ");
817 printf("Expand an existing container ... ");
818 }
819 signals->Expand();
820 }
793ff80c 821 }
e23fbb27 822 // The same for the dictionary
0c24ba98 823 if (!fSimpleSim) {
824 for (iDict = 0; iDict < kNDict; iDict++) {
825 dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
826 if (dictionary[iDict]->GetNtime() == 0) {
827 dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
828 }
829 else {
830 if (fCompress) dictionary[iDict]->Expand();
831 }
e23fbb27 832 }
833 }
834 if (fDebug > 1) printf("done\n");
835 detectorOld = detector;
793ff80c 836 }
837
e23fbb27 838 // Rotate the sectors on top of each other
0c24ba98 839 if (fSimpleSim) {
840 rot[0] = pos[0];
841 rot[1] = pos[1];
842 rot[2] = pos[2];
843 }
844 else {
845 fGeo->Rotate(detector,pos,rot);
846 }
e23fbb27 847
848 // The driftlength. It is negative if the hit is in the
849 // amplification region.
850 Float_t driftlength = time0 - rot[0];
851
852 // Take also the drift in the amplification region into account
853 // The drift length is at the moment still the same, regardless of
854 // the position relativ to the wire. This non-isochronity needs still
855 // to be implemented.
856 Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
4487dad0 857 if (fPar->ExBOn()) driftlengthL /= TMath::Sqrt(fPar->GetLorentzFactor());
e23fbb27 858
859 // Loop over all electrons of this hit
860 // TR photons produce hits with negative charge
861 Int_t nEl = ((Int_t) TMath::Abs(q));
862 for (Int_t iEl = 0; iEl < nEl; iEl++) {
863
864 xyz[0] = rot[0];
865 xyz[1] = rot[1];
866 xyz[2] = rot[2];
867
868 // Electron attachment
4487dad0 869 if (fPar->ElAttachOn()) {
e23fbb27 870 if (gRandom->Rndm() < (driftlengthL * elAttachProp))
871 continue;
872 }
873
874 // Apply the diffusion smearing
4487dad0 875 if (fPar->DiffusionOn()) {
876 if (!(fPar->Diffusion(driftlengthL,xyz))) continue;
e23fbb27 877 }
f7336fa3 878
e23fbb27 879 // Apply E x B effects (depends on drift direction)
4487dad0 880 if (fPar->ExBOn()) {
881 if (!(fPar->ExB(driftlength+kAmWidth,xyz))) continue;
e23fbb27 882 }
793ff80c 883
e23fbb27 884 // The electron position after diffusion and ExB in pad coordinates
885 // The pad row (z-direction)
886 Float_t rowDist = xyz[2] - row0;
887 Int_t rowE = ((Int_t) (rowDist * divideRow));
888 if ((rowE < 0) || (rowE >= nRowMax)) continue;
889 Float_t rowOffset = ((((Float_t) rowE) + 0.5) * rowPadSize) - rowDist;
890
891 // The pad column (rphi-direction)
4487dad0 892 Float_t col0tilt = fPar->Col0Tilted(col0,rowOffset,plane);
e23fbb27 893 Float_t colDist = xyz[1] - col0tilt;
894 Int_t colE = ((Int_t) (colDist * divideCol));
895 if ((colE < 0) || (colE >= nColMax)) continue;
896 Float_t colOffset = ((((Float_t) colE) + 0.5) * colPadSize) - colDist;
897
898 // The time bin (negative for hits in the amplification region)
899 // In the amplification region the electrons drift from both sides
900 // to the middle (anode wire plane)
901 Float_t timeDist = time0 - xyz[0];
902 Float_t timeOffset = 0;
903 Int_t timeE = 0;
904 if (timeDist > 0) {
905 // The time bin
906 timeE = ((Int_t) (timeDist * divideTime));
907 // The distance of the position to the middle of the timebin
908 timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
909 }
910 else {
911 // Difference between half of the amplification gap width and
912 // the distance to the anode wire
913 Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
914 // The time bin
915 timeE = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
916 // The distance of the position to the middle of the timebin
917 timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
918 }
872a7aba 919
e23fbb27 920 // Apply the gas gain including fluctuations
921 Float_t ggRndm = 0.0;
922 do {
923 ggRndm = gRandom->Rndm();
924 } while (ggRndm <= 0);
4487dad0 925 Int_t signal = (Int_t) (-fPar->GetGasGain() * TMath::Log(ggRndm));
e23fbb27 926
927 // Apply the pad response
4487dad0 928 if (fPar->PRFOn()) {
e23fbb27 929 // The distance of the electron to the center of the pad
930 // in units of pad width
931 Float_t dist = - colOffset * divideCol;
4487dad0 932 if (!(fPar->PadResponse(signal,dist,plane,padSignal))) continue;
793ff80c 933 }
e23fbb27 934 else {
935 padSignal[0] = 0.0;
936 padSignal[1] = signal;
937 padSignal[2] = 0.0;
938 }
939
940 // Sample the time response inside the drift region
941 // + additional time bins before and after.
942 // The sampling is done always in the middle of the time bin
943 for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg, -nTimeBefore)
944 ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter )
945 ;iTimeBin++) {
946
947 // Apply the time response
948 Float_t timeResponse = 1.0;
949 Float_t crossTalk = 0.0;
950 Float_t time = (iTimeBin - timeE) * timeBinSize + timeOffset;
4487dad0 951 if (fPar->TRFOn()) {
952 timeResponse = fPar->TimeResponse(time);
e23fbb27 953 }
4487dad0 954 if (fPar->CTOn()) {
955 crossTalk = fPar->CrossTalk(time);
e23fbb27 956 }
957
958 signalOld[0] = 0.0;
959 signalOld[1] = 0.0;
960 signalOld[2] = 0.0;
f7336fa3 961
e23fbb27 962 for (iPad = 0; iPad < kNpad; iPad++) {
963
964 Int_t colPos = colE + iPad - 1;
965 if (colPos < 0) continue;
966 if (colPos >= nColMax) break;
967
968 // Add the signals
969 // Note: The time bin number is shifted by nTimeBefore to avoid negative
970 // time bins. This has to be subtracted later.
971 Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
972 signalOld[iPad] = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
973 if( colPos != colE ) {
974 signalOld[iPad] += padSignal[iPad] * (timeResponse + crossTalk);
975 }
976 else {
977 signalOld[iPad] += padSignal[iPad] * timeResponse;
978 }
979 signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
980
981 // Store the track index in the dictionary
982 // Note: We store index+1 in order to allow the array to be compressed
0c24ba98 983 if ((signalOld[iPad] > 0) && (!fSimpleSim)) {
e23fbb27 984 for (iDict = 0; iDict < kNDict; iDict++) {
985 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
986 ,colPos
987 ,iCurrentTimeBin);
988 if (oldTrack == track+1) break;
989 if (oldTrack == 0) {
990 dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
991 break;
992 }
71d9fa7b 993 }
994 }
872a7aba 995
e23fbb27 996 } // Loop: pads
f7336fa3 997
e23fbb27 998 } // Loop: time bins
f7336fa3 999
e23fbb27 1000 } // Loop: electrons of a single hit
f7336fa3 1001
e23fbb27 1002 } // If: detector and test hit
f7336fa3 1003
e23fbb27 1004 hit = (AliTRDhit *) fTRD->NextHit();
f7336fa3 1005
e23fbb27 1006 } // Loop: hits of one primary track
1007
1008 } // Loop: primary tracks
1009
1010 if (fDebug > 0) {
1011 printf("<AliTRDdigitizer::MakeDigits> ");
abaf1f1d 1012 printf("Finished analyzing %d hits\n",countHits);
1013 }
793ff80c 1014
0c24ba98 1015 // The coupling factor
1016 Float_t coupling = fPar->GetPadCoupling()
1017 * fPar->GetTimeCoupling();
1018
1019 // The conversion factor
1020 Float_t convert = kEl2fC
1021 * fPar->GetChipGain();
6244debe 1022
793ff80c 1023 // Loop through all chambers to finalize the digits
0c24ba98 1024 Int_t iDetBeg = 0;
1025 Int_t iDetEnd = AliTRDgeometry::Ndet();
1026 if (fSimpleSim) {
1027 iDetBeg = fSimpleDet;
1028 iDetEnd = iDetBeg + 1;
1029 }
1030 for (Int_t iDet = iDetBeg; iDet < iDetEnd; iDet++) {
793ff80c 1031
872a7aba 1032 Int_t plane = fGeo->GetPlane(iDet);
1033 Int_t sector = fGeo->GetSector(iDet);
1034 Int_t chamber = fGeo->GetChamber(iDet);
4487dad0 1035 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
1036 Int_t nColMax = fPar->GetColMax(plane);
1037 Int_t nTimeMax = fPar->GetTimeMax();
1038 Int_t nTimeTotal = fPar->GetTimeTotal();
793ff80c 1039
e23fbb27 1040 Double_t *inADC = new Double_t[nTimeTotal];
1041 Double_t *outADC = new Double_t[nTimeTotal];
1042
1043 if (fDebug > 0) {
1044 printf("<AliTRDdigitizer::MakeDigits> ");
793ff80c 1045 printf("Digitization for chamber %d\n",iDet);
1046 }
da581aea 1047
793ff80c 1048 // Add a container for the digits of this detector
abaf1f1d 1049 digits = fDigitsManager->GetDigits(iDet);
793ff80c 1050 // Allocate memory space for the digits buffer
0c24ba98 1051 if (digits->GetNtime() == 0) {
1052 digits->Allocate(nRowMax,nColMax,nTimeTotal);
1053 }
1054 else if (fSimpleSim) {
1055 digits->Clear();
1056 }
1057
793ff80c 1058 // Get the signal container
1059 signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1060 if (signals->GetNtime() == 0) {
1061 // Create missing containers
872a7aba 1062 signals->Allocate(nRowMax,nColMax,nTimeTotal);
793ff80c 1063 }
1064 else {
1065 // Expand the container if neccessary
1066 if (fCompress) signals->Expand();
1067 }
1068 // Create the missing dictionary containers
0c24ba98 1069 if (!fSimpleSim) {
1070 for (iDict = 0; iDict < kNDict; iDict++) {
1071 dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1072 if (dictionary[iDict]->GetNtime() == 0) {
1073 dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1074 }
1075 }
793ff80c 1076 }
f7336fa3 1077
793ff80c 1078 Int_t nDigits = 0;
1079
6244debe 1080 // Don't create noise in detectors that are switched off
1081 if (CheckDetector(plane,chamber,sector)) {
1082
1083 // Create the digits for this chamber
872a7aba 1084 for (iRow = 0; iRow < nRowMax; iRow++ ) {
1085 for (iCol = 0; iCol < nColMax; iCol++ ) {
6244debe 1086
e23fbb27 1087 // Create summable digits
1088 if (fSDigits) {
6244debe 1089
e23fbb27 1090 for (iTime = 0; iTime < nTimeTotal; iTime++) {
872a7aba 1091 Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
abaf1f1d 1092 signalAmp *= fSDigitsScale;
e23fbb27 1093 signalAmp = TMath::Min(signalAmp,(Float_t) 1.0e9);
abaf1f1d 1094 Int_t adc = (Int_t) signalAmp;
4487dad0 1095 if (adc > 0) nDigits++;
872a7aba 1096 digits->SetDataUnchecked(iRow,iCol,iTime,adc);
c1e4b257 1097 }
6244debe 1098
e23fbb27 1099 }
1100 // Create normal digits
1101 else {
6244debe 1102
e23fbb27 1103 for (iTime = 0; iTime < nTimeTotal; iTime++) {
1104 Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
0c24ba98 1105 // Pad and time coupling
1106 signalAmp *= coupling;
6244debe 1107 // Add the noise
4487dad0 1108 signalAmp = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fPar->GetNoise()),0.0);
6244debe 1109 // Convert to mV
1110 signalAmp *= convert;
1111 // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1112 // signal is larger than fADCinRange
1113 Int_t adc = 0;
4487dad0 1114 if (signalAmp >= fPar->GetADCinRange()) {
1115 adc = ((Int_t) fPar->GetADCoutRange());
6244debe 1116 }
1117 else {
4487dad0 1118 adc = ((Int_t) (signalAmp * (fPar->GetADCoutRange()
1119 / fPar->GetADCinRange())));
6244debe 1120 }
e23fbb27 1121 inADC[iTime] = adc;
1122 outADC[iTime] = adc;
1123 }
6244debe 1124
e23fbb27 1125 // Apply the tail cancelation via the digital filter
4487dad0 1126 if (fPar->TCOn()) {
1127 DeConvExp(inADC,outADC,nTimeTotal,fPar->GetTCnexp());
e23fbb27 1128 }
1129
1130 for (iTime = 0; iTime < nTimeTotal; iTime++) {
6244debe 1131 // Store the amplitude of the digit if above threshold
4487dad0 1132 if (outADC[iTime] > fPar->GetADCthreshold()) {
e23fbb27 1133 if (fDebug > 2) {
1134 printf(" iRow = %d, iCol = %d, iTime = %d, adc = %f\n"
1135 ,iRow,iCol,iTime,outADC[iTime]);
6244debe 1136 }
1137 nDigits++;
3583ac1a 1138 digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
6244debe 1139 }
6244debe 1140 }
1141
1142 }
e23fbb27 1143
6244debe 1144 }
793ff80c 1145 }
6244debe 1146
793ff80c 1147 }
1148
1149 // Compress the arrays
0c24ba98 1150 if (!fSimpleSim) {
1151 digits->Compress(1,0);
1152 for (iDict = 0; iDict < kNDict; iDict++) {
1153 dictionary[iDict]->Compress(1,0);
1154 }
f7336fa3 1155
0c24ba98 1156 totalSizeDigits += digits->GetSize();
1157 totalSizeDict0 += dictionary[0]->GetSize();
1158 totalSizeDict1 += dictionary[1]->GetSize();
1159 totalSizeDict2 += dictionary[2]->GetSize();
f7336fa3 1160
0c24ba98 1161 Float_t nPixel = nRowMax * nColMax * nTimeMax;
1162 if (fDebug > 0) {
1163 printf("<AliTRDdigitizer::MakeDigits> ");
1164 printf("Found %d digits in detector %d (%3.0f).\n"
1165 ,nDigits,iDet
1166 ,100.0 * ((Float_t) nDigits) / nPixel);
1167 }
1168
1169 if (fCompress) signals->Compress(1,0);
abaf1f1d 1170
0c24ba98 1171 }
f7336fa3 1172
e23fbb27 1173 delete [] inADC;
1174 delete [] outADC;
1175
f7336fa3 1176 }
1177
0c24ba98 1178 if (signalsArray) {
1179 delete signalsArray;
1180 signalsArray = 0;
1181 }
1182
e23fbb27 1183 if (fDebug > 0) {
1184 printf("<AliTRDdigitizer::MakeDigits> ");
abaf1f1d 1185 printf("Total number of analyzed hits = %d\n",countHits);
0c24ba98 1186 if (!fSimpleSim) {
1187 printf("<AliTRDdigitizer::MakeDigits> ");
1188 printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
1189 ,totalSizeDict0
1190 ,totalSizeDict1
1191 ,totalSizeDict2);
1192 }
abaf1f1d 1193 }
1194
1195 return kTRUE;
1196
1197}
1198
1199//_____________________________________________________________________________
1200void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1201{
1202 //
1203 // Add a digits manager for s-digits to the input list.
1204 //
1205
1206 fSDigitsManagerList->Add(man);
1207
1208}
1209
1210//_____________________________________________________________________________
b1113c6b 1211void AliTRDdigitizer::DeleteSDigitsManager()
1212{
1213 //
1214 // Removes digits manager from the input list.
1215 //
1216
1217 fSDigitsManagerList->Delete();
1218
1219}
1220
1221//_____________________________________________________________________________
abaf1f1d 1222Bool_t AliTRDdigitizer::ConvertSDigits()
1223{
1224 //
1225 // Converts s-digits to normal digits
1226 //
1227
1228 // Number of track dictionary arrays
1229 const Int_t kNDict = AliTRDdigitsManager::kNDict;
1230
1231 // Converts number of electrons to fC
1232 const Double_t kEl2fC = 1.602E-19 * 1.0E15;
1233
1234 Int_t iDict = 0;
e23fbb27 1235 Int_t iRow;
1236 Int_t iCol;
1237 Int_t iTime;
abaf1f1d 1238
4487dad0 1239 if (!fPar) {
1240 fPar = new AliTRDparameter("TRDparameter","Standard parameter");
1241 if (fDebug > 0) {
1242 printf("<AliTRDdigitizer::ConvertSDigits> ");
1243 printf("Create the default parameter object\n");
1244 }
abaf1f1d 1245 }
1246
1247 Double_t sDigitsScale = 1.0 / GetSDigitsScale();
4487dad0 1248 Double_t noise = fPar->GetNoise();
1249 Double_t padCoupling = fPar->GetPadCoupling();
1250 Double_t timeCoupling = fPar->GetTimeCoupling();
1251 Double_t chipGain = fPar->GetChipGain();
0c24ba98 1252 Double_t coupling = padCoupling * timeCoupling;
1253 Double_t convert = kEl2fC * chipGain;
4487dad0 1254 Double_t adcInRange = fPar->GetADCinRange();
1255 Double_t adcOutRange = fPar->GetADCoutRange();
1256 Int_t adcThreshold = fPar->GetADCthreshold();
abaf1f1d 1257
1258 AliTRDdataArrayI *digitsIn;
1259 AliTRDdataArrayI *digitsOut;
1260 AliTRDdataArrayI *dictionaryIn[kNDict];
1261 AliTRDdataArrayI *dictionaryOut[kNDict];
1262
1263 // Loop through the detectors
1264 for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
da581aea 1265
e23fbb27 1266 if (fDebug > 0) {
1267 printf("<AliTRDdigitizer::ConvertSDigits> ");
abaf1f1d 1268 printf("Convert detector %d to digits.\n",iDet);
1269 }
1270
1271 Int_t plane = fGeo->GetPlane(iDet);
1272 Int_t sector = fGeo->GetSector(iDet);
1273 Int_t chamber = fGeo->GetChamber(iDet);
4487dad0 1274 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
1275 Int_t nColMax = fPar->GetColMax(plane);
1276 Int_t nTimeTotal = fPar->GetTimeTotal();
abaf1f1d 1277
e23fbb27 1278 Double_t *inADC = new Double_t[nTimeTotal];
1279 Double_t *outADC = new Double_t[nTimeTotal];
1280
abaf1f1d 1281 digitsIn = fSDigitsManager->GetDigits(iDet);
1282 digitsIn->Expand();
1283 digitsOut = fDigitsManager->GetDigits(iDet);
1284 digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1285 for (iDict = 0; iDict < kNDict; iDict++) {
1286 dictionaryIn[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1287 dictionaryIn[iDict]->Expand();
1288 dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1289 dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1290 }
1291
e23fbb27 1292 for (iRow = 0; iRow < nRowMax; iRow++ ) {
1293 for (iCol = 0; iCol < nColMax; iCol++ ) {
abaf1f1d 1294
e23fbb27 1295 for (iTime = 0; iTime < nTimeTotal; iTime++) {
abaf1f1d 1296 Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1297 signal *= sDigitsScale;
0c24ba98 1298 // Pad and time coupling
1299 signal *= coupling;
abaf1f1d 1300 // Add the noise
1301 signal = TMath::Max((Double_t) gRandom->Gaus(signal,noise),0.0);
1302 // Convert to mV
1303 signal *= convert;
1304 // Convert to ADC counts. Set the overflow-bit adcOutRange if the
1305 // signal is larger than adcInRange
1306 Int_t adc = 0;
1307 if (signal >= adcInRange) {
1308 adc = ((Int_t) adcOutRange);
1309 }
1310 else {
1311 adc = ((Int_t) (signal * (adcOutRange / adcInRange)));
1312 }
e23fbb27 1313 inADC[iTime] = adc;
1314 outADC[iTime] = adc;
1315 }
1316
1317 // Apply the tail cancelation via the digital filter
4487dad0 1318 if (fPar->TCOn()) {
1319 DeConvExp(inADC,outADC,nTimeTotal,fPar->GetTCnexp());
e23fbb27 1320 }
1321
1322 for (iTime = 0; iTime < nTimeTotal; iTime++) {
abaf1f1d 1323 // Store the amplitude of the digit if above threshold
e23fbb27 1324 if (outADC[iTime] > adcThreshold) {
3583ac1a 1325 digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
e23fbb27 1326 // Copy the dictionary
1327 for (iDict = 0; iDict < kNDict; iDict++) {
1328 Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1329 dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
1330 }
abaf1f1d 1331 }
abaf1f1d 1332 }
e23fbb27 1333
abaf1f1d 1334 }
1335 }
1336
1337 if (fCompress) {
1338 digitsIn->Compress(1,0);
1339 digitsOut->Compress(1,0);
1340 for (iDict = 0; iDict < kNDict; iDict++) {
1341 dictionaryIn[iDict]->Compress(1,0);
1342 dictionaryOut[iDict]->Compress(1,0);
1343 }
1344 }
1345
e23fbb27 1346 delete [] inADC;
1347 delete [] outADC;
1348
abaf1f1d 1349 }
f7336fa3 1350
1351 return kTRUE;
1352
1353}
1354
1355//_____________________________________________________________________________
abaf1f1d 1356Bool_t AliTRDdigitizer::MergeSDigits()
16bf9884 1357{
1358 //
abaf1f1d 1359 // Merges the input s-digits:
1360 // - The amplitude of the different inputs are summed up.
1361 // - Of the track IDs from the input dictionaries only one is
1362 // kept for each input. This works for maximal 3 different merged inputs.
16bf9884 1363 //
1364
abaf1f1d 1365 // Number of track dictionary arrays
1366 const Int_t kNDict = AliTRDdigitsManager::kNDict;
1367
4487dad0 1368 if (!fPar) {
1369 fPar = new AliTRDparameter("TRDparameter","Standard parameter");
1370 if (fDebug > 0) {
1371 printf("<AliTRDdigitizer::MergeSDigits> ");
1372 printf("Create the default parameter object\n");
1373 }
1374 }
1375
abaf1f1d 1376 Int_t iDict = 0;
e23fbb27 1377 Int_t jDict = 0;
abaf1f1d 1378
1379 AliTRDdataArrayI *digitsA;
1380 AliTRDdataArrayI *digitsB;
1381 AliTRDdataArrayI *dictionaryA[kNDict];
1382 AliTRDdataArrayI *dictionaryB[kNDict];
1383
1384 // Get the first s-digits
1385 fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1386 if (!fSDigitsManager) return kFALSE;
1387
1388 // Loop through the other sets of s-digits
1389 AliTRDdigitsManager *mergeSDigitsManager;
1390 mergeSDigitsManager = (AliTRDdigitsManager *)
1391 fSDigitsManagerList->After(fSDigitsManager);
1392
e23fbb27 1393 if (fDebug > 0) {
abaf1f1d 1394 if (mergeSDigitsManager) {
e23fbb27 1395 printf("<AliTRDdigitizer::MergeSDigits> ");
1396 printf("Merge %d input files.\n",fSDigitsManagerList->GetSize());
abaf1f1d 1397 }
1398 else {
e23fbb27 1399 printf("<AliTRDdigitizer::MergeSDigits> ");
abaf1f1d 1400 printf("Only one input file.\n");
1401 }
1402 }
1403
1404 Int_t iMerge = 0;
1405 while (mergeSDigitsManager) {
1406
1407 iMerge++;
1408
1409 // Loop through the detectors
1410 for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1411
1412 Int_t plane = fGeo->GetPlane(iDet);
1413 Int_t sector = fGeo->GetSector(iDet);
1414 Int_t chamber = fGeo->GetChamber(iDet);
4487dad0 1415 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
1416 Int_t nColMax = fPar->GetColMax(plane);
1417 Int_t nTimeTotal = fPar->GetTimeTotal();
abaf1f1d 1418
1419 // Loop through the pixels of one detector and add the signals
1420 digitsA = fSDigitsManager->GetDigits(iDet);
1421 digitsB = mergeSDigitsManager->GetDigits(iDet);
1422 digitsA->Expand();
1423 digitsB->Expand();
1424 for (iDict = 0; iDict < kNDict; iDict++) {
1425 dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1426 dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
1427 dictionaryA[iDict]->Expand();
1428 dictionaryB[iDict]->Expand();
1429 }
1430
4487dad0 1431 // Merge only detectors that contain a signal
1432 Bool_t doMerge = kTRUE;
1433 if (fMergeSignalOnly) {
1434 if (digitsA->GetOverThreshold(0) == 0) {
1435 doMerge = kFALSE;
1436 }
abaf1f1d 1437 }
1438
4487dad0 1439 if (doMerge) {
abaf1f1d 1440
4487dad0 1441 if (fDebug > 0) {
1442 printf("<AliTRDdigitizer::MergeSDigits> ");
1443 printf("Merge detector %d of input no.%d\n",iDet,iMerge+1);
1444 }
abaf1f1d 1445
4487dad0 1446 for (Int_t iRow = 0; iRow < nRowMax; iRow++ ) {
1447 for (Int_t iCol = 0; iCol < nColMax; iCol++ ) {
1448 for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1449
1450 // Add the amplitudes of the summable digits
1451 Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
1452 Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
1453 ampA += ampB;
1454 digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
1455
1456 // Add the mask to the track id if defined.
1457 for (iDict = 0; iDict < kNDict; iDict++) {
1458 Int_t trackB = dictionaryB[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1459 if ((fMasks) && (trackB > 0)) {
1460 for (jDict = 0; jDict < kNDict; jDict++) {
1461 Int_t trackA = dictionaryA[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1462 if (trackA == 0) {
1463 trackA = trackB + fMasks[iMerge];
1464 dictionaryA[iDict]->SetDataUnchecked(iRow,iCol,iTime,trackA);
1465 }
1466 }
1467 }
e23fbb27 1468 }
abaf1f1d 1469
4487dad0 1470 }
abaf1f1d 1471 }
4487dad0 1472 }
1473
abaf1f1d 1474 }
1475
1476 if (fCompress) {
1477 digitsA->Compress(1,0);
1478 digitsB->Compress(1,0);
1479 for (iDict = 0; iDict < kNDict; iDict++) {
1480 dictionaryA[iDict]->Compress(1,0);
1481 dictionaryB[iDict]->Compress(1,0);
1482 }
1483 }
1484
1485 }
1486
1487 // The next set of s-digits
1488 mergeSDigitsManager = (AliTRDdigitsManager *)
1489 fSDigitsManagerList->After(mergeSDigitsManager);
1490
1491 }
1492
16bf9884 1493 return kTRUE;
1494
1495}
1496
1497//_____________________________________________________________________________
abaf1f1d 1498Bool_t AliTRDdigitizer::SDigits2Digits()
1499{
1500 //
1501 // Merges the input s-digits and converts them to normal digits
1502 //
1503
1504 if (!MergeSDigits()) return kFALSE;
1505
1506 return ConvertSDigits();
1507
1508}
1509
1510//_____________________________________________________________________________
793ff80c 1511Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
1512{
1513 //
1514 // Checks whether a detector is enabled
1515 //
1516
0c24ba98 1517 if (fSimpleSim) return kTRUE;
1518
793ff80c 1519 if ((fTRD->GetSensChamber() >= 0) &&
1520 (fTRD->GetSensChamber() != chamber)) return kFALSE;
1521 if ((fTRD->GetSensPlane() >= 0) &&
c1e4b257 1522 (fTRD->GetSensPlane() != plane)) return kFALSE;
793ff80c 1523 if ( fTRD->GetSensSector() >= 0) {
1524 Int_t sens1 = fTRD->GetSensSector();
1525 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
1526 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
1527 * AliTRDgeometry::Nsect();
1528 if (sens1 < sens2) {
1529 if ((sector < sens1) || (sector >= sens2)) return kFALSE;
1530 }
1531 else {
1532 if ((sector < sens1) && (sector >= sens2)) return kFALSE;
1533 }
1534 }
1535
1536 return kTRUE;
1537
1538}
1539
1540//_____________________________________________________________________________
0a29d0f1 1541Bool_t AliTRDdigitizer::WriteDigits() const
f7336fa3 1542{
1543 //
1544 // Writes out the TRD-digits and the dictionaries
1545 //
1546
da581aea 1547 // Store the digits and the dictionary in the tree
abaf1f1d 1548 return fDigitsManager->WriteDigits();
f7336fa3 1549
1550}
793ff80c 1551
1552//_____________________________________________________________________________
e23fbb27 1553void AliTRDdigitizer::DeConvExp(Double_t *source, Double_t *target
1554 , Int_t n, Int_t nexp)
1555{
1556 //
1557 // Does the deconvolution by the digital filter.
1558 //
1559 // Author: Marcus Gutfleisch, KIP Heidelberg
1560 // Optimized for: New TRF from Venelin Angelov, simulated with CADENCE
1561 // Pad-ground capacitance = 25 pF
1562 // Pad-pad cross talk capacitance = 6 pF
1563 // For 10 MHz digitization, corresponding to 20 time bins
1564 // in the drift region
1565 //
1566
1567 Double_t rates[2];
1568 Double_t coefficients[2];
1569
1570 /* initialize (coefficient = alpha, rates = lambda) */
1571
1572 if( nexp == 1 ) {
1573 rates[0] = 0.466998;
1574 /* no rescaling */
1575 coefficients[0] = 1.0;
1576 }
1577 if( nexp == 2 ) {
1578 rates[0] = 0.8988162;
1579 coefficients[0] = 0.11392069;
1580 rates[1] = 0.3745688;
1581 coefficients[1] = 0.8860793;
1582 /* no rescaling */
1583 Float_t sumc = coefficients[0]+coefficients[1];
1584 coefficients[0] /= sumc;
1585 coefficients[1] /= sumc;
1586 }
1587
1588 Int_t i, k;
1589 Double_t reminder[2];
1590 Double_t correction, result;
1591
1592 /* attention: computation order is important */
1593 correction=0.0;
1594 for ( k = 0; k < nexp; k++ ) reminder[k]=0.0;
1595
1596 for ( i = 0; i < n; i++ ) {
1597 result = ( source[i] - correction ); /* no rescaling */
1598 target[i] = result;
1599
1600 for ( k = 0; k < nexp; k++ ) reminder[k] = rates[k]
1601 * ( reminder[k] + coefficients[k] * result);
1602
1603 correction=0.0;
1604 for ( k = 0; k < nexp; k++ ) correction += reminder[k];
1605 }
1606
1607}
1608
8e64dd77 1609//_____________________________________________________________________________
1610void AliTRDdigitizer::InitOutput(TFile *file, Int_t iEvent)
1611{
1612 //
1613 // Initializes the output branches
1614 //
1615
1616 fEvent = iEvent;
1617 fDigitsManager->MakeTreeAndBranches(file,iEvent);
1618
1619}