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