1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.37 2002/04/29 11:50:47 cblume
19 Change initialization of gAlice in the merging case
21 Revision 1.36 2002/04/12 12:13:23 cblume
24 Revision 1.35 2002/03/28 14:59:07 cblume
27 Revision 1.34 2002/03/25 20:00:44 cblume
28 Introduce parameter class and regions of interest for merging
30 Revision 1.33 2002/02/12 17:32:03 cblume
31 Rearrange the deleting of the list of sdigitsmanager
33 Revision 1.32 2002/02/12 16:07:21 cblume
36 Revision 1.31 2002/02/11 14:27:11 cblume
37 New pad plane design, new TRF+PRF, tail cancelation, cross talk
39 Revision 1.30 2001/11/19 08:44:08 cblume
40 Fix bugs reported by Rene
42 Revision 1.29 2001/11/14 19:44:25 hristov
43 Numeric const casted (Alpha)
45 Revision 1.28 2001/11/14 16:35:58 cblume
46 Inherits now from AliDetector
48 Revision 1.27 2001/11/14 10:50:45 cblume
49 Changes in digits IO. Add merging of summable digits
51 Revision 1.26 2001/11/06 17:19:41 cblume
52 Add detailed geometry and simple simulator
54 Revision 1.25 2001/06/27 09:54:44 cblume
55 Moved fField initialization to InitDetector()
57 Revision 1.24 2001/05/21 16:45:47 hristov
58 Last minute changes (C.Blume)
60 Revision 1.23 2001/05/07 08:04:48 cblume
61 New TRF and PRF. Speedup of the code. Digits from amplification region included
63 Revision 1.22 2001/03/30 14:40:14 cblume
64 Update of the digitization parameter
66 Revision 1.21 2001/03/13 09:30:35 cblume
67 Update of digitization. Moved digit branch definition to AliTRD
69 Revision 1.20 2001/02/25 20:19:00 hristov
70 Minor correction: loop variable declared only once for HP, Sun
72 Revision 1.19 2001/02/14 18:22:26 cblume
73 Change in the geometry of the padplane
75 Revision 1.18 2001/01/26 19:56:57 hristov
76 Major upgrade of AliRoot code
78 Revision 1.17 2000/12/08 12:53:27 cblume
79 Change in Copy() function for HP-compiler
81 Revision 1.16 2000/12/07 12:20:46 cblume
82 Go back to array compression. Use sampled PRF to speed up digitization
84 Revision 1.15 2000/11/23 14:34:08 cblume
85 Fixed bug in expansion routine of arrays (initialize buffers properly)
87 Revision 1.14 2000/11/20 08:54:44 cblume
88 Switch off compression as default
90 Revision 1.13 2000/11/10 14:57:52 cblume
91 Changes in the geometry constants for the DEC compiler
93 Revision 1.12 2000/11/01 14:53:20 cblume
94 Merge with TRD-develop
96 Revision 1.1.4.9 2000/10/26 17:00:22 cblume
97 Fixed bug in CheckDetector()
99 Revision 1.1.4.8 2000/10/23 13:41:35 cblume
100 Added protection against Log(0) in the gas gain calulation
102 Revision 1.1.4.7 2000/10/17 02:27:34 cblume
103 Get rid of global constants
105 Revision 1.1.4.6 2000/10/16 01:16:53 cblume
106 Changed timebin 0 to be the one closest to the readout
108 Revision 1.1.4.5 2000/10/15 23:34:29 cblume
109 Faster version of the digitizer
111 Revision 1.1.4.4 2000/10/06 16:49:46 cblume
114 Revision 1.1.4.3 2000/10/04 16:34:58 cblume
115 Replace include files by forward declarations
117 Revision 1.1.4.2 2000/09/22 14:41:10 cblume
118 Bug fix in PRF. Included time response. New structure
120 Revision 1.10 2000/10/05 07:27:53 cblume
121 Changes in the header-files by FCA
123 Revision 1.9 2000/10/02 21:28:19 fca
124 Removal of useless dependecies via forward declarations
126 Revision 1.8 2000/06/09 11:10:07 cblume
127 Compiler warnings and coding conventions, next round
129 Revision 1.7 2000/06/08 18:32:58 cblume
130 Make code compliant to coding conventions
132 Revision 1.6 2000/06/07 16:27:32 cblume
133 Try to remove compiler warnings on Sun and HP
135 Revision 1.5 2000/05/09 16:38:57 cblume
136 Removed PadResponse(). Merge problem
138 Revision 1.4 2000/05/08 15:53:45 cblume
139 Resolved merge conflict
141 Revision 1.3 2000/04/28 14:49:27 cblume
142 Only one declaration of iDict in MakeDigits()
144 Revision 1.1.4.1 2000/05/08 14:42:04 cblume
145 Introduced AliTRDdigitsManager
147 Revision 1.1 2000/02/28 19:00:13 cblume
152 ///////////////////////////////////////////////////////////////////////////////
154 // Creates and handles digits from TRD hits //
155 // Author: C. Blume (C.Blume@gsi.de) //
157 // The following effects are included: //
160 // - Gas gain including fluctuations //
161 // - Pad-response (simple Gaussian approximation) //
162 // - Time-response //
163 // - Electronics noise //
164 // - Electronics gain //
166 // - ADC threshold //
167 // The corresponding parameter can be adjusted via the various //
168 // Set-functions. If these parameters are not explicitly set, default //
169 // values are used (see Init-function). //
170 // As an example on how to use this class to produce digits from hits //
171 // have a look at the macro hits2digits.C //
172 // The production of summable digits is demonstrated in hits2sdigits.C //
173 // and the subsequent conversion of the s-digits into normal digits is //
174 // explained in sdigits2digits.C. //
176 ///////////////////////////////////////////////////////////////////////////////
192 #include "AliRunDigitizer.h"
195 #include "AliTRDhit.h"
196 #include "AliTRDdigitizer.h"
197 #include "AliTRDdataArrayI.h"
198 #include "AliTRDdataArrayF.h"
199 #include "AliTRDsegmentArray.h"
200 #include "AliTRDdigitsManager.h"
201 #include "AliTRDgeometry.h"
202 #include "AliTRDparameter.h"
204 ClassImp(AliTRDdigitizer)
206 //_____________________________________________________________________________
207 AliTRDdigitizer::AliTRDdigitizer()
210 // AliTRDdigitizer default constructor
215 fSDigitsManagerList = 0;
226 fMergeSignalOnly = kFALSE;
230 //_____________________________________________________________________________
231 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)
232 :AliDigitizer(name,title)
235 // AliTRDdigitizer constructor
240 fSDigitsManagerList = 0;
250 fMergeSignalOnly = kFALSE;
252 // For the summable digits
253 fSDigitsScale = 100.;
257 //_____________________________________________________________________________
258 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
259 , const Text_t *name, const Text_t *title)
260 :AliDigitizer(manager,name,title)
263 // AliTRDdigitizer constructor
268 fSDigitsManagerList = 0;
278 fMergeSignalOnly = kFALSE;
280 // For the summable digits
281 fSDigitsScale = 100.;
285 //_____________________________________________________________________________
286 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
287 :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
290 // AliTRDdigitizer constructor
295 fSDigitsManagerList = 0;
305 fMergeSignalOnly = kFALSE;
307 // For the summable digits
308 fSDigitsScale = 100.;
312 //_____________________________________________________________________________
313 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
316 // AliTRDdigitizer copy constructor
319 ((AliTRDdigitizer &) d).Copy(*this);
323 //_____________________________________________________________________________
324 AliTRDdigitizer::~AliTRDdigitizer()
327 // AliTRDdigitizer destructor
336 if (fDigitsManager) {
337 delete fDigitsManager;
341 if (fSDigitsManager) {
342 delete fSDigitsManager;
346 if (fSDigitsManagerList) {
347 delete fSDigitsManagerList;
348 fSDigitsManagerList = 0;
358 //_____________________________________________________________________________
359 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
362 // Assignment operator
365 if (this != &d) ((AliTRDdigitizer &) d).Copy(*this);
370 //_____________________________________________________________________________
371 void AliTRDdigitizer::Copy(TObject &d)
377 ((AliTRDdigitizer &) d).fInputFile = 0;
378 ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
379 ((AliTRDdigitizer &) d).fSDigitsManager = 0;
380 ((AliTRDdigitizer &) d).fDigitsManager = 0;
381 ((AliTRDdigitizer &) d).fTRD = 0;
382 ((AliTRDdigitizer &) d).fGeo = 0;
383 ((AliTRDdigitizer &) d).fMasks = 0;
384 ((AliTRDdigitizer &) d).fEvent = 0;
385 ((AliTRDdigitizer &) d).fPar = 0;
386 ((AliTRDdigitizer &) d).fCompress = fCompress;
387 ((AliTRDdigitizer &) d).fDebug = fDebug ;
388 ((AliTRDdigitizer &) d).fSDigits = fSDigits;
389 ((AliTRDdigitizer &) d).fSDigitsScale = fSDigitsScale;
390 ((AliTRDdigitizer &) d).fMergeSignalOnly = fMergeSignalOnly;
394 //_____________________________________________________________________________
395 void AliTRDdigitizer::Exec(Option_t* option)
398 // Executes the merging
403 AliTRDdigitsManager *sdigitsManager;
405 TString optionString = option;
406 if (optionString.Contains("deb")) {
408 if (optionString.Contains("2")) {
411 printf("<AliTRDdigitizer::Exec> ");
412 printf("Called with debug option %d\n",fDebug);
415 // The AliRoot file is already connected by the manager
418 printf("<AliTRDdigitizer::Exec> ");
419 printf("AliRun object found on file.\n");
423 printf("<AliTRDdigitizer::Exec> ");
424 printf("Could not find AliRun object.\n");
428 Int_t nInput = fManager->GetNinputs();
429 fMasks = new Int_t[nInput];
430 for (iInput = 0; iInput < nInput; iInput++) {
431 fMasks[iInput] = fManager->GetMask(iInput);
437 for (iInput = 0; iInput < nInput; iInput++) {
440 printf("<AliTRDdigitizer::Exec> ");
441 printf("Add input stream %d\n",iInput);
444 // check if the input tree exists
445 if (!fManager->GetInputTreeTRDS(iInput)) {
446 printf("<AliTRDdigitizer::Exec> ");
447 printf("Input stream %d does not exist\n",iInput);
451 // Read the s-digits via digits manager
452 sdigitsManager = new AliTRDdigitsManager();
453 sdigitsManager->SetDebug(fDebug);
454 sdigitsManager->SetSDigits(kTRUE);
455 sdigitsManager->ReadDigits(fManager->GetInputTreeTRDS(iInput));
457 // Add the s-digits to the input list
458 AddSDigitsManager(sdigitsManager);
462 // Convert the s-digits to normal digits
464 printf("<AliTRDdigitizer::Exec> ");
465 printf("Do the conversion\n");
471 printf("<AliTRDdigitizer::Exec> ");
472 printf("Write the digits\n");
474 fDigitsManager->MakeBranch(fManager->GetTreeDTRD());
475 fDigitsManager->WriteDigits();
477 printf("<AliTRDdigitizer::Exec> ");
481 DeleteSDigitsManager();
485 //_____________________________________________________________________________
486 Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
489 // Opens a ROOT-file with TRD-hits and reads in the hit-tree
492 // Connect the AliRoot file containing Geometry, Kine, and Hits
493 fInputFile = (TFile *) gROOT->GetListOfFiles()->FindObject(file);
496 printf("<AliTRDdigitizer::Open> ");
497 printf("Open the AliROOT-file %s.\n",file);
499 fInputFile = new TFile(file,"UPDATE");
503 printf("<AliTRDdigitizer::Open> ");
504 printf("%s is already open.\n",file);
508 gAlice = (AliRun *) fInputFile->Get("gAlice");
511 printf("<AliTRDdigitizer::Open> ");
512 printf("AliRun object found on file.\n");
516 printf("<AliTRDdigitizer::Open> ");
517 printf("Could not find AliRun object.\n");
523 // Import the Trees for the event nEvent in the file
524 Int_t nparticles = gAlice->GetEvent(fEvent);
525 if (nparticles <= 0) {
526 printf("<AliTRDdigitizer::Open> ");
527 printf("No entries in the trees for event %d.\n",fEvent);
531 if (InitDetector()) {
540 //_____________________________________________________________________________
541 Bool_t AliTRDdigitizer::InitDetector()
544 // Sets the pointer to the TRD detector and the geometry
547 // Get the pointer to the detector class and check for version 1
548 fTRD = (AliTRD *) gAlice->GetDetector("TRD");
550 printf("<AliTRDdigitizer::InitDetector> ");
551 printf("No TRD module found\n");
554 if (fTRD->IsVersion() != 1) {
555 printf("<AliTRDdigitizer::InitDetector> ");
556 printf("TRD must be version 1 (slow simulator).\n");
561 fGeo = fTRD->GetGeometry();
563 printf("<AliTRDdigitizer::InitDetector> ");
564 printf("Geometry version %d\n",fGeo->IsVersion());
567 // Create a digits manager
568 fDigitsManager = new AliTRDdigitsManager();
569 fDigitsManager->SetSDigits(fSDigits);
570 fDigitsManager->CreateArrays();
571 fDigitsManager->SetEvent(fEvent);
572 fDigitsManager->SetDebug(fDebug);
574 // The list for the input s-digits manager to be merged
575 fSDigitsManagerList = new TList();
581 //_____________________________________________________________________________
582 Bool_t AliTRDdigitizer::MakeBranch(const Char_t *file) const
585 // Create the branches for the digits array
588 return fDigitsManager->MakeBranch(file);
592 //_____________________________________________________________________________
593 Bool_t AliTRDdigitizer::MakeDigits()
599 ///////////////////////////////////////////////////////////////
601 ///////////////////////////////////////////////////////////////
603 // Converts number of electrons to fC
604 const Double_t kEl2fC = 1.602E-19 * 1.0E15;
606 ///////////////////////////////////////////////////////////////
608 // Number of pads included in the pad response
609 const Int_t kNpad = 3;
611 // Number of track dictionary arrays
612 const Int_t kNDict = AliTRDdigitsManager::kNDict;
614 // Half the width of the amplification region
615 const Float_t kAmWidth = AliTRDgeometry::AmThick() / 2.;
617 Int_t iRow, iCol, iTime, iPad;
621 Int_t totalSizeDigits = 0;
622 Int_t totalSizeDict0 = 0;
623 Int_t totalSizeDict1 = 0;
624 Int_t totalSizeDict2 = 0;
626 Int_t timeTRDbeg = 0;
627 Int_t timeTRDend = 1;
632 Float_t padSignal[kNpad];
633 Float_t signalOld[kNpad];
635 AliTRDdataArrayF *signals = 0;
636 AliTRDdataArrayI *digits = 0;
637 AliTRDdataArrayI *dictionary[kNDict];
639 // Create a default parameter class if none is defined
641 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
643 printf("<AliTRDdigitizer::MakeDigits> ");
644 printf("Create the default parameter object\n");
648 // Create a container for the amplitudes
649 AliTRDsegmentArray *signalsArray
650 = new AliTRDsegmentArray("AliTRDdataArrayF"
651 ,AliTRDgeometry::Ndet());
654 timeTRDbeg = ((Int_t) (-fPar->GetTRFlo() / fPar->GetTimeBinSize())) - 1;
655 timeTRDend = ((Int_t) ( fPar->GetTRFhi() / fPar->GetTimeBinSize())) - 1;
657 printf("<AliTRDdigitizer::MakeDigits> ");
658 printf("Sample the TRF between -%d and %d\n",timeTRDbeg,timeTRDend);
662 Float_t elAttachProp = fPar->GetElAttachProp() / 100.;
665 printf("<AliTRDdigitizer::MakeDigits> ");
666 printf("No geometry defined\n");
671 printf("<AliTRDdigitizer::MakeDigits> ");
672 printf("Start creating digits.\n");
675 // Get the pointer to the hit tree
676 TTree *hitTree = gAlice->TreeH();
678 // Get the number of entries in the hit tree
679 // (Number of primary particles creating a hit somewhere)
680 Int_t nTrack = (Int_t) hitTree->GetEntries();
682 printf("<AliTRDdigitizer::MakeDigits> ");
683 printf("Found %d primary particles\n",nTrack);
686 Int_t detectorOld = -1;
689 // Loop through all entries in the tree
690 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
693 nBytes += hitTree->GetEvent(iTrack);
695 // Loop through the TRD hits
697 AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
706 Float_t q = hit->GetCharge();
707 Int_t track = hit->Track();
708 Int_t detector = hit->GetDetector();
709 Int_t plane = fGeo->GetPlane(detector);
710 Int_t sector = fGeo->GetSector(detector);
711 Int_t chamber = fGeo->GetChamber(detector);
712 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
713 Int_t nColMax = fPar->GetColMax(plane);
714 Int_t nTimeMax = fPar->GetTimeMax();
715 Int_t nTimeBefore = fPar->GetTimeBefore();
716 Int_t nTimeAfter = fPar->GetTimeAfter();
717 Int_t nTimeTotal = fPar->GetTimeTotal();
718 Float_t row0 = fPar->GetRow0(plane,chamber,sector);
719 Float_t col0 = fPar->GetCol0(plane);
720 Float_t time0 = fPar->GetTime0(plane);
721 Float_t rowPadSize = fPar->GetRowPadSize(plane,chamber,sector);
722 Float_t colPadSize = fPar->GetColPadSize(plane);
723 Float_t timeBinSize = fPar->GetTimeBinSize();
724 Float_t divideRow = 1.0 / rowPadSize;
725 Float_t divideCol = 1.0 / colPadSize;
726 Float_t divideTime = 1.0 / timeBinSize;
729 printf("Analyze hit no. %d ",iHit);
730 printf("-----------------------------------------------------------\n");
732 printf("plane = %d, sector = %d, chamber = %d\n"
733 ,plane,sector,chamber);
734 printf("nRowMax = %d, nColMax = %d, nTimeMax = %d\n"
735 ,nRowMax,nColMax,nTimeMax);
736 printf("nTimeBefore = %d, nTimeAfter = %d, nTimeTotal = %d\n"
737 ,nTimeBefore,nTimeAfter,nTimeTotal);
738 printf("row0 = %f, col0 = %f, time0 = %f\n"
740 printf("rowPadSize = %f, colPadSize = %f, timeBinSize = %f\n"
741 ,rowPadSize,colPadSize,timeBinSize);
744 // Don't analyze test hits and switched off detectors
745 if ((CheckDetector(plane,chamber,sector)) &&
746 (((Int_t) q) != 0)) {
748 if (detector != detectorOld) {
751 printf("<AliTRDdigitizer::MakeDigits> ");
752 printf("Get new container. New det = %d, Old det = %d\n"
753 ,detector,detectorOld);
755 // Compress the old one if enabled
756 if ((fCompress) && (detectorOld > -1)) {
758 printf("<AliTRDdigitizer::MakeDigits> ");
759 printf("Compress the old container ...");
761 signals->Compress(1,0);
762 for (iDict = 0; iDict < kNDict; iDict++) {
763 dictionary[iDict]->Compress(1,0);
765 if (fDebug > 1) printf("done\n");
767 // Get the new container
768 signals = (AliTRDdataArrayF *) signalsArray->At(detector);
769 if (signals->GetNtime() == 0) {
770 // Allocate a new one if not yet existing
772 printf("<AliTRDdigitizer::MakeDigits> ");
773 printf("Allocate a new container ... ");
775 signals->Allocate(nRowMax,nColMax,nTimeTotal);
778 // Expand an existing one
781 printf("<AliTRDdigitizer::MakeDigits> ");
782 printf("Expand an existing container ... ");
787 // The same for the dictionary
788 for (iDict = 0; iDict < kNDict; iDict++) {
789 dictionary[iDict] = fDigitsManager->GetDictionary(detector,iDict);
790 if (dictionary[iDict]->GetNtime() == 0) {
791 dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
794 if (fCompress) dictionary[iDict]->Expand();
797 if (fDebug > 1) printf("done\n");
798 detectorOld = detector;
801 // Rotate the sectors on top of each other
802 fGeo->Rotate(detector,pos,rot);
804 // The driftlength. It is negative if the hit is in the
805 // amplification region.
806 Float_t driftlength = time0 - rot[0];
808 // Take also the drift in the amplification region into account
809 // The drift length is at the moment still the same, regardless of
810 // the position relativ to the wire. This non-isochronity needs still
811 // to be implemented.
812 Float_t driftlengthL = TMath::Abs(driftlength + kAmWidth);
813 if (fPar->ExBOn()) driftlengthL /= TMath::Sqrt(fPar->GetLorentzFactor());
815 // Loop over all electrons of this hit
816 // TR photons produce hits with negative charge
817 Int_t nEl = ((Int_t) TMath::Abs(q));
818 for (Int_t iEl = 0; iEl < nEl; iEl++) {
824 // Electron attachment
825 if (fPar->ElAttachOn()) {
826 if (gRandom->Rndm() < (driftlengthL * elAttachProp))
830 // Apply the diffusion smearing
831 if (fPar->DiffusionOn()) {
832 if (!(fPar->Diffusion(driftlengthL,xyz))) continue;
835 // Apply E x B effects (depends on drift direction)
837 if (!(fPar->ExB(driftlength+kAmWidth,xyz))) continue;
840 // The electron position after diffusion and ExB in pad coordinates
841 // The pad row (z-direction)
842 Float_t rowDist = xyz[2] - row0;
843 Int_t rowE = ((Int_t) (rowDist * divideRow));
844 if ((rowE < 0) || (rowE >= nRowMax)) continue;
845 Float_t rowOffset = ((((Float_t) rowE) + 0.5) * rowPadSize) - rowDist;
847 // The pad column (rphi-direction)
848 Float_t col0tilt = fPar->Col0Tilted(col0,rowOffset,plane);
849 Float_t colDist = xyz[1] - col0tilt;
850 Int_t colE = ((Int_t) (colDist * divideCol));
851 if ((colE < 0) || (colE >= nColMax)) continue;
852 Float_t colOffset = ((((Float_t) colE) + 0.5) * colPadSize) - colDist;
854 // The time bin (negative for hits in the amplification region)
855 // In the amplification region the electrons drift from both sides
856 // to the middle (anode wire plane)
857 Float_t timeDist = time0 - xyz[0];
858 Float_t timeOffset = 0;
862 timeE = ((Int_t) (timeDist * divideTime));
863 // The distance of the position to the middle of the timebin
864 timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) - timeDist;
867 // Difference between half of the amplification gap width and
868 // the distance to the anode wire
869 Float_t anodeDist = kAmWidth - TMath::Abs(timeDist + kAmWidth);
871 timeE = -1 * (((Int_t ) (anodeDist * divideTime)) + 1);
872 // The distance of the position to the middle of the timebin
873 timeOffset = ((((Float_t) timeE) + 0.5) * timeBinSize) + anodeDist;
876 // Apply the gas gain including fluctuations
877 Float_t ggRndm = 0.0;
879 ggRndm = gRandom->Rndm();
880 } while (ggRndm <= 0);
881 Int_t signal = (Int_t) (-fPar->GetGasGain() * TMath::Log(ggRndm));
883 // Apply the pad response
885 // The distance of the electron to the center of the pad
886 // in units of pad width
887 Float_t dist = - colOffset * divideCol;
888 if (!(fPar->PadResponse(signal,dist,plane,padSignal))) continue;
892 padSignal[1] = signal;
896 // Sample the time response inside the drift region
897 // + additional time bins before and after.
898 // The sampling is done always in the middle of the time bin
899 for (Int_t iTimeBin = TMath::Max(timeE-timeTRDbeg, -nTimeBefore)
900 ;iTimeBin < TMath::Min(timeE+timeTRDend,nTimeMax+nTimeAfter )
903 // Apply the time response
904 Float_t timeResponse = 1.0;
905 Float_t crossTalk = 0.0;
906 Float_t time = (iTimeBin - timeE) * timeBinSize + timeOffset;
908 timeResponse = fPar->TimeResponse(time);
911 crossTalk = fPar->CrossTalk(time);
918 for (iPad = 0; iPad < kNpad; iPad++) {
920 Int_t colPos = colE + iPad - 1;
921 if (colPos < 0) continue;
922 if (colPos >= nColMax) break;
925 // Note: The time bin number is shifted by nTimeBefore to avoid negative
926 // time bins. This has to be subtracted later.
927 Int_t iCurrentTimeBin = iTimeBin + nTimeBefore;
928 signalOld[iPad] = signals->GetDataUnchecked(rowE,colPos,iCurrentTimeBin);
929 if( colPos != colE ) {
930 signalOld[iPad] += padSignal[iPad] * (timeResponse + crossTalk);
933 signalOld[iPad] += padSignal[iPad] * timeResponse;
935 signals->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,signalOld[iPad]);
937 // Store the track index in the dictionary
938 // Note: We store index+1 in order to allow the array to be compressed
939 if (signalOld[iPad] > 0) {
940 for (iDict = 0; iDict < kNDict; iDict++) {
941 Int_t oldTrack = dictionary[iDict]->GetDataUnchecked(rowE
944 if (oldTrack == track+1) break;
946 dictionary[iDict]->SetDataUnchecked(rowE,colPos,iCurrentTimeBin,track+1);
956 } // Loop: electrons of a single hit
958 } // If: detector and test hit
960 hit = (AliTRDhit *) fTRD->NextHit();
962 } // Loop: hits of one primary track
964 } // Loop: primary tracks
967 printf("<AliTRDdigitizer::MakeDigits> ");
968 printf("Finished analyzing %d hits\n",countHits);
971 // The total conversion factor
972 Float_t convert = kEl2fC * fPar->GetPadCoupling()
973 * fPar->GetTimeCoupling()
974 * fPar->GetChipGain();
976 // Loop through all chambers to finalize the digits
977 for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
979 Int_t plane = fGeo->GetPlane(iDet);
980 Int_t sector = fGeo->GetSector(iDet);
981 Int_t chamber = fGeo->GetChamber(iDet);
982 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
983 Int_t nColMax = fPar->GetColMax(plane);
984 Int_t nTimeMax = fPar->GetTimeMax();
985 Int_t nTimeTotal = fPar->GetTimeTotal();
987 Double_t *inADC = new Double_t[nTimeTotal];
988 Double_t *outADC = new Double_t[nTimeTotal];
991 printf("<AliTRDdigitizer::MakeDigits> ");
992 printf("Digitization for chamber %d\n",iDet);
995 // Add a container for the digits of this detector
996 digits = fDigitsManager->GetDigits(iDet);
997 // Allocate memory space for the digits buffer
998 digits->Allocate(nRowMax,nColMax,nTimeTotal);
1000 // Get the signal container
1001 signals = (AliTRDdataArrayF *) signalsArray->At(iDet);
1002 if (signals->GetNtime() == 0) {
1003 // Create missing containers
1004 signals->Allocate(nRowMax,nColMax,nTimeTotal);
1007 // Expand the container if neccessary
1008 if (fCompress) signals->Expand();
1010 // Create the missing dictionary containers
1011 for (iDict = 0; iDict < kNDict; iDict++) {
1012 dictionary[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1013 if (dictionary[iDict]->GetNtime() == 0) {
1014 dictionary[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1020 // Don't create noise in detectors that are switched off
1021 if (CheckDetector(plane,chamber,sector)) {
1023 // Create the digits for this chamber
1024 for (iRow = 0; iRow < nRowMax; iRow++ ) {
1025 for (iCol = 0; iCol < nColMax; iCol++ ) {
1027 // Create summable digits
1030 for (iTime = 0; iTime < nTimeTotal; iTime++) {
1031 Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1032 signalAmp *= fSDigitsScale;
1033 signalAmp = TMath::Min(signalAmp,(Float_t) 1.0e9);
1034 Int_t adc = (Int_t) signalAmp;
1035 if (adc > 0) nDigits++;
1036 digits->SetDataUnchecked(iRow,iCol,iTime,adc);
1040 // Create normal digits
1043 for (iTime = 0; iTime < nTimeTotal; iTime++) {
1044 Float_t signalAmp = signals->GetDataUnchecked(iRow,iCol,iTime);
1046 signalAmp = TMath::Max((Double_t) gRandom->Gaus(signalAmp,fPar->GetNoise()),0.0);
1048 signalAmp *= convert;
1049 // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1050 // signal is larger than fADCinRange
1052 if (signalAmp >= fPar->GetADCinRange()) {
1053 adc = ((Int_t) fPar->GetADCoutRange());
1056 adc = ((Int_t) (signalAmp * (fPar->GetADCoutRange()
1057 / fPar->GetADCinRange())));
1060 outADC[iTime] = adc;
1063 // Apply the tail cancelation via the digital filter
1065 DeConvExp(inADC,outADC,nTimeTotal,fPar->GetTCnexp());
1068 for (iTime = 0; iTime < nTimeTotal; iTime++) {
1069 // Store the amplitude of the digit if above threshold
1070 if (outADC[iTime] > fPar->GetADCthreshold()) {
1072 printf(" iRow = %d, iCol = %d, iTime = %d, adc = %f\n"
1073 ,iRow,iCol,iTime,outADC[iTime]);
1076 digits->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1087 // Compress the arrays
1088 digits->Compress(1,0);
1089 for (iDict = 0; iDict < kNDict; iDict++) {
1090 dictionary[iDict]->Compress(1,0);
1093 totalSizeDigits += digits->GetSize();
1094 totalSizeDict0 += dictionary[0]->GetSize();
1095 totalSizeDict1 += dictionary[1]->GetSize();
1096 totalSizeDict2 += dictionary[2]->GetSize();
1098 Float_t nPixel = nRowMax * nColMax * nTimeMax;
1100 printf("<AliTRDdigitizer::MakeDigits> ");
1101 printf("Found %d digits in detector %d (%3.0f).\n"
1103 ,100.0 * ((Float_t) nDigits) / nPixel);
1106 if (fCompress) signals->Compress(1,0);
1114 printf("<AliTRDdigitizer::MakeDigits> ");
1115 printf("Total number of analyzed hits = %d\n",countHits);
1116 printf("<AliTRDdigitizer::MakeDigits> ");
1117 printf("Total digits data size = %d, %d, %d, %d\n",totalSizeDigits
1127 //_____________________________________________________________________________
1128 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
1131 // Add a digits manager for s-digits to the input list.
1134 fSDigitsManagerList->Add(man);
1138 //_____________________________________________________________________________
1139 void AliTRDdigitizer::DeleteSDigitsManager()
1142 // Removes digits manager from the input list.
1145 fSDigitsManagerList->Delete();
1149 //_____________________________________________________________________________
1150 Bool_t AliTRDdigitizer::ConvertSDigits()
1153 // Converts s-digits to normal digits
1156 // Number of track dictionary arrays
1157 const Int_t kNDict = AliTRDdigitsManager::kNDict;
1159 // Converts number of electrons to fC
1160 const Double_t kEl2fC = 1.602E-19 * 1.0E15;
1168 fPar = new AliTRDparameter("TRDparameter","Standard parameter");
1170 printf("<AliTRDdigitizer::ConvertSDigits> ");
1171 printf("Create the default parameter object\n");
1175 Double_t sDigitsScale = 1.0 / GetSDigitsScale();
1176 Double_t noise = fPar->GetNoise();
1177 Double_t padCoupling = fPar->GetPadCoupling();
1178 Double_t timeCoupling = fPar->GetTimeCoupling();
1179 Double_t chipGain = fPar->GetChipGain();
1180 Double_t convert = kEl2fC * padCoupling * timeCoupling * chipGain;;
1181 Double_t adcInRange = fPar->GetADCinRange();
1182 Double_t adcOutRange = fPar->GetADCoutRange();
1183 Int_t adcThreshold = fPar->GetADCthreshold();
1185 AliTRDdataArrayI *digitsIn;
1186 AliTRDdataArrayI *digitsOut;
1187 AliTRDdataArrayI *dictionaryIn[kNDict];
1188 AliTRDdataArrayI *dictionaryOut[kNDict];
1190 // Loop through the detectors
1191 for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1194 printf("<AliTRDdigitizer::ConvertSDigits> ");
1195 printf("Convert detector %d to digits.\n",iDet);
1198 Int_t plane = fGeo->GetPlane(iDet);
1199 Int_t sector = fGeo->GetSector(iDet);
1200 Int_t chamber = fGeo->GetChamber(iDet);
1201 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
1202 Int_t nColMax = fPar->GetColMax(plane);
1203 Int_t nTimeTotal = fPar->GetTimeTotal();
1205 Double_t *inADC = new Double_t[nTimeTotal];
1206 Double_t *outADC = new Double_t[nTimeTotal];
1208 digitsIn = fSDigitsManager->GetDigits(iDet);
1210 digitsOut = fDigitsManager->GetDigits(iDet);
1211 digitsOut->Allocate(nRowMax,nColMax,nTimeTotal);
1212 for (iDict = 0; iDict < kNDict; iDict++) {
1213 dictionaryIn[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1214 dictionaryIn[iDict]->Expand();
1215 dictionaryOut[iDict] = fDigitsManager->GetDictionary(iDet,iDict);
1216 dictionaryOut[iDict]->Allocate(nRowMax,nColMax,nTimeTotal);
1219 for (iRow = 0; iRow < nRowMax; iRow++ ) {
1220 for (iCol = 0; iCol < nColMax; iCol++ ) {
1222 for (iTime = 0; iTime < nTimeTotal; iTime++) {
1223 Double_t signal = (Double_t) digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1224 signal *= sDigitsScale;
1226 signal = TMath::Max((Double_t) gRandom->Gaus(signal,noise),0.0);
1229 // Convert to ADC counts. Set the overflow-bit adcOutRange if the
1230 // signal is larger than adcInRange
1232 if (signal >= adcInRange) {
1233 adc = ((Int_t) adcOutRange);
1236 adc = ((Int_t) (signal * (adcOutRange / adcInRange)));
1239 outADC[iTime] = adc;
1242 // Apply the tail cancelation via the digital filter
1244 DeConvExp(inADC,outADC,nTimeTotal,fPar->GetTCnexp());
1247 for (iTime = 0; iTime < nTimeTotal; iTime++) {
1248 // Store the amplitude of the digit if above threshold
1249 if (outADC[iTime] > adcThreshold) {
1250 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1251 // Copy the dictionary
1252 for (iDict = 0; iDict < kNDict; iDict++) {
1253 Int_t track = dictionaryIn[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1254 dictionaryOut[iDict]->SetDataUnchecked(iRow,iCol,iTime,track);
1263 digitsIn->Compress(1,0);
1264 digitsOut->Compress(1,0);
1265 for (iDict = 0; iDict < kNDict; iDict++) {
1266 dictionaryIn[iDict]->Compress(1,0);
1267 dictionaryOut[iDict]->Compress(1,0);
1280 //_____________________________________________________________________________
1281 Bool_t AliTRDdigitizer::MergeSDigits()
1284 // Merges the input s-digits:
1285 // - The amplitude of the different inputs are summed up.
1286 // - Of the track IDs from the input dictionaries only one is
1287 // kept for each input. This works for maximal 3 different merged inputs.
1290 // Number of track dictionary arrays
1291 const Int_t kNDict = AliTRDdigitsManager::kNDict;
1294 fPar = new AliTRDparameter("TRDparameter","Standard parameter");
1296 printf("<AliTRDdigitizer::MergeSDigits> ");
1297 printf("Create the default parameter object\n");
1304 AliTRDdataArrayI *digitsA;
1305 AliTRDdataArrayI *digitsB;
1306 AliTRDdataArrayI *dictionaryA[kNDict];
1307 AliTRDdataArrayI *dictionaryB[kNDict];
1309 // Get the first s-digits
1310 fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1311 if (!fSDigitsManager) return kFALSE;
1313 // Loop through the other sets of s-digits
1314 AliTRDdigitsManager *mergeSDigitsManager;
1315 mergeSDigitsManager = (AliTRDdigitsManager *)
1316 fSDigitsManagerList->After(fSDigitsManager);
1319 if (mergeSDigitsManager) {
1320 printf("<AliTRDdigitizer::MergeSDigits> ");
1321 printf("Merge %d input files.\n",fSDigitsManagerList->GetSize());
1324 printf("<AliTRDdigitizer::MergeSDigits> ");
1325 printf("Only one input file.\n");
1330 while (mergeSDigitsManager) {
1334 // Loop through the detectors
1335 for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1337 Int_t plane = fGeo->GetPlane(iDet);
1338 Int_t sector = fGeo->GetSector(iDet);
1339 Int_t chamber = fGeo->GetChamber(iDet);
1340 Int_t nRowMax = fPar->GetRowMax(plane,chamber,sector);
1341 Int_t nColMax = fPar->GetColMax(plane);
1342 Int_t nTimeTotal = fPar->GetTimeTotal();
1344 // Loop through the pixels of one detector and add the signals
1345 digitsA = fSDigitsManager->GetDigits(iDet);
1346 digitsB = mergeSDigitsManager->GetDigits(iDet);
1349 for (iDict = 0; iDict < kNDict; iDict++) {
1350 dictionaryA[iDict] = fSDigitsManager->GetDictionary(iDet,iDict);
1351 dictionaryB[iDict] = mergeSDigitsManager->GetDictionary(iDet,iDict);
1352 dictionaryA[iDict]->Expand();
1353 dictionaryB[iDict]->Expand();
1356 // Merge only detectors that contain a signal
1357 Bool_t doMerge = kTRUE;
1358 if (fMergeSignalOnly) {
1359 if (digitsA->GetOverThreshold(0) == 0) {
1367 printf("<AliTRDdigitizer::MergeSDigits> ");
1368 printf("Merge detector %d of input no.%d\n",iDet,iMerge+1);
1371 for (Int_t iRow = 0; iRow < nRowMax; iRow++ ) {
1372 for (Int_t iCol = 0; iCol < nColMax; iCol++ ) {
1373 for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1375 // Add the amplitudes of the summable digits
1376 Int_t ampA = digitsA->GetDataUnchecked(iRow,iCol,iTime);
1377 Int_t ampB = digitsB->GetDataUnchecked(iRow,iCol,iTime);
1379 digitsA->SetDataUnchecked(iRow,iCol,iTime,ampA);
1381 // Add the mask to the track id if defined.
1382 for (iDict = 0; iDict < kNDict; iDict++) {
1383 Int_t trackB = dictionaryB[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1384 if ((fMasks) && (trackB > 0)) {
1385 for (jDict = 0; jDict < kNDict; jDict++) {
1386 Int_t trackA = dictionaryA[iDict]->GetDataUnchecked(iRow,iCol,iTime);
1388 trackA = trackB + fMasks[iMerge];
1389 dictionaryA[iDict]->SetDataUnchecked(iRow,iCol,iTime,trackA);
1402 digitsA->Compress(1,0);
1403 digitsB->Compress(1,0);
1404 for (iDict = 0; iDict < kNDict; iDict++) {
1405 dictionaryA[iDict]->Compress(1,0);
1406 dictionaryB[iDict]->Compress(1,0);
1412 // The next set of s-digits
1413 mergeSDigitsManager = (AliTRDdigitsManager *)
1414 fSDigitsManagerList->After(mergeSDigitsManager);
1422 //_____________________________________________________________________________
1423 Bool_t AliTRDdigitizer::SDigits2Digits()
1426 // Merges the input s-digits and converts them to normal digits
1429 if (!MergeSDigits()) return kFALSE;
1431 return ConvertSDigits();
1435 //_____________________________________________________________________________
1436 Bool_t AliTRDdigitizer::CheckDetector(Int_t plane, Int_t chamber, Int_t sector)
1439 // Checks whether a detector is enabled
1442 if ((fTRD->GetSensChamber() >= 0) &&
1443 (fTRD->GetSensChamber() != chamber)) return kFALSE;
1444 if ((fTRD->GetSensPlane() >= 0) &&
1445 (fTRD->GetSensPlane() != plane)) return kFALSE;
1446 if ( fTRD->GetSensSector() >= 0) {
1447 Int_t sens1 = fTRD->GetSensSector();
1448 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
1449 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
1450 * AliTRDgeometry::Nsect();
1451 if (sens1 < sens2) {
1452 if ((sector < sens1) || (sector >= sens2)) return kFALSE;
1455 if ((sector < sens1) && (sector >= sens2)) return kFALSE;
1463 //_____________________________________________________________________________
1464 Bool_t AliTRDdigitizer::WriteDigits() const
1467 // Writes out the TRD-digits and the dictionaries
1470 // Store the digits and the dictionary in the tree
1471 return fDigitsManager->WriteDigits();
1475 //_____________________________________________________________________________
1476 void AliTRDdigitizer::DeConvExp(Double_t *source, Double_t *target
1477 , Int_t n, Int_t nexp)
1480 // Does the deconvolution by the digital filter.
1482 // Author: Marcus Gutfleisch, KIP Heidelberg
1483 // Optimized for: New TRF from Venelin Angelov, simulated with CADENCE
1484 // Pad-ground capacitance = 25 pF
1485 // Pad-pad cross talk capacitance = 6 pF
1486 // For 10 MHz digitization, corresponding to 20 time bins
1487 // in the drift region
1491 Double_t coefficients[2];
1493 /* initialize (coefficient = alpha, rates = lambda) */
1496 rates[0] = 0.466998;
1498 coefficients[0] = 1.0;
1501 rates[0] = 0.8988162;
1502 coefficients[0] = 0.11392069;
1503 rates[1] = 0.3745688;
1504 coefficients[1] = 0.8860793;
1506 Float_t sumc = coefficients[0]+coefficients[1];
1507 coefficients[0] /= sumc;
1508 coefficients[1] /= sumc;
1512 Double_t reminder[2];
1513 Double_t correction, result;
1515 /* attention: computation order is important */
1517 for ( k = 0; k < nexp; k++ ) reminder[k]=0.0;
1519 for ( i = 0; i < n; i++ ) {
1520 result = ( source[i] - correction ); /* no rescaling */
1523 for ( k = 0; k < nexp; k++ ) reminder[k] = rates[k]
1524 * ( reminder[k] + coefficients[k] * result);
1527 for ( k = 0; k < nexp; k++ ) correction += reminder[k];
1532 //_____________________________________________________________________________
1533 void AliTRDdigitizer::InitOutput(TFile *file, Int_t iEvent)
1536 // Initializes the output branches
1540 fDigitsManager->MakeTreeAndBranches(file,iEvent);