003e6b292c7358002644460996060aaab7465fe2
[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 //                                                                        //
22 //  Authors: C. Blume (blume@ikf.uni-frankfurt.de)                        //
23 //           C. Lippmann                                                  //
24 //           B. Vulpescu                                                  //
25 //                                                                        //
26 //  The following effects are included:                                   //
27 //      - Diffusion                                                       //
28 //      - ExB effects                                                     //
29 //      - Gas gain including fluctuations                                 //
30 //      - Pad-response (simple Gaussian approximation)                    //
31 //      - Time-response                                                   //
32 //      - Electronics noise                                               //
33 //      - Electronics gain                                                //
34 //      - Digitization                                                    //
35 //      - Zero suppression                                                //
36 //                                                                        //
37 ////////////////////////////////////////////////////////////////////////////
38
39 #include <TGeoManager.h>
40 #include <TList.h>
41 #include <TMath.h>
42 #include <TRandom.h>
43 #include <TTree.h>
44
45 #include "AliRun.h"
46 #include "AliMC.h"
47 #include "AliRunLoader.h"
48 #include "AliLoader.h"
49 #include "AliConfig.h"
50 #include "AliRunDigitizer.h"
51 #include "AliRunLoader.h"
52 #include "AliLoader.h"
53 #include "AliLog.h"
54
55 #include "AliTRD.h"
56 #include "AliTRDhit.h"
57 #include "AliTRDdigitizer.h"
58 #include "AliTRDarrayDictionary.h"
59 #include "AliTRDarrayADC.h"
60 #include "AliTRDarraySignal.h"
61 #include "AliTRDdigitsManager.h"
62 #include "AliTRDgeometry.h"
63 #include "AliTRDpadPlane.h"
64 #include "AliTRDcalibDB.h"
65 #include "AliTRDSimParam.h"
66 #include "AliTRDCommonParam.h"
67 #include "AliTRDfeeParam.h"
68 #include "AliTRDmcmSim.h"
69 #include "Cal/AliTRDCalROC.h"
70 #include "Cal/AliTRDCalDet.h"
71
72 ClassImp(AliTRDdigitizer)
73
74 //_____________________________________________________________________________
75 AliTRDdigitizer::AliTRDdigitizer()
76   :AliDigitizer()
77   ,fRunLoader(0)
78   ,fDigitsManager(0)
79   ,fSDigitsManager(0)
80   ,fSDigitsManagerList(0)
81   ,fTRD(0)
82   ,fGeo(0)
83   ,fEvent(0)
84   ,fMasks(0)
85   ,fCompress(kTRUE)
86   ,fSDigits(kFALSE)
87   ,fMergeSignalOnly(kFALSE)
88 {
89   //
90   // AliTRDdigitizer default constructor
91   //
92   
93   Init();
94
95 }
96
97 //_____________________________________________________________________________
98 AliTRDdigitizer::AliTRDdigitizer(const Text_t *name, const Text_t *title)              
99   :AliDigitizer(name,title)
100   ,fRunLoader(0)
101   ,fDigitsManager(0)
102   ,fSDigitsManager(0)
103   ,fSDigitsManagerList(0)
104   ,fTRD(0)
105   ,fGeo(0)
106   ,fEvent(0)
107   ,fMasks(0)
108   ,fCompress(kTRUE)
109   ,fSDigits(kFALSE)
110   ,fMergeSignalOnly(kFALSE)
111 {
112   //
113   // AliTRDdigitizer constructor
114   //
115
116   Init();
117
118 }
119
120 //_____________________________________________________________________________
121 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager
122                                , const Text_t *name, const Text_t *title)
123   :AliDigitizer(manager,name,title)
124   ,fRunLoader(0)
125   ,fDigitsManager(0)
126   ,fSDigitsManager(0)
127   ,fSDigitsManagerList(0)
128   ,fTRD(0)
129   ,fGeo(0)
130   ,fEvent(0)
131   ,fMasks(0)
132   ,fCompress(kTRUE)
133   ,fSDigits(kFALSE)
134   ,fMergeSignalOnly(kFALSE)
135 {
136   //
137   // AliTRDdigitizer constructor
138   //
139
140   Init();
141
142 }
143
144 //_____________________________________________________________________________
145 AliTRDdigitizer::AliTRDdigitizer(AliRunDigitizer *manager)
146   :AliDigitizer(manager,"AliTRDdigitizer","TRD digitizer")
147   ,fRunLoader(0)
148   ,fDigitsManager(0)
149   ,fSDigitsManager(0)
150   ,fSDigitsManagerList(0)
151   ,fTRD(0)
152   ,fGeo(0)
153   ,fEvent(0)
154   ,fMasks(0)
155   ,fCompress(kTRUE)
156   ,fSDigits(kFALSE)
157   ,fMergeSignalOnly(kFALSE)
158 {
159   //
160   // AliTRDdigitizer constructor
161   //
162
163   Init();
164
165 }
166
167 //_____________________________________________________________________________
168 Bool_t AliTRDdigitizer::Init()
169 {
170   //
171   // Initialize the digitizer with default values
172   //
173
174   fRunLoader          = 0;
175   fDigitsManager      = 0;
176   fSDigitsManager     = 0;
177   fSDigitsManagerList = 0;
178   fTRD                = 0;
179   fGeo                = 0;
180
181   fEvent              = 0;
182   fMasks              = 0;
183   fCompress           = kTRUE;
184   fSDigits            = kFALSE;
185   fMergeSignalOnly    = kFALSE;
186   
187   return AliDigitizer::Init();
188
189 }
190
191 //_____________________________________________________________________________
192 AliTRDdigitizer::AliTRDdigitizer(const AliTRDdigitizer &d)
193   :AliDigitizer(d)
194   ,fRunLoader(0)
195   ,fDigitsManager(0)
196   ,fSDigitsManager(0)
197   ,fSDigitsManagerList(0)
198   ,fTRD(0)
199   ,fGeo(0)
200   ,fEvent(0)
201   ,fMasks(0)
202   ,fCompress(d.fCompress)
203   ,fSDigits(d.fSDigits)
204   ,fMergeSignalOnly(d.fMergeSignalOnly)
205 {
206   //
207   // AliTRDdigitizer copy constructor
208   //
209
210 }
211
212 //_____________________________________________________________________________
213 AliTRDdigitizer::~AliTRDdigitizer()
214 {
215   //
216   // AliTRDdigitizer destructor
217   //
218
219   if (fDigitsManager) {
220     delete fDigitsManager;
221     fDigitsManager      = 0;
222   }
223
224   if (fDigitsManager) {              //typo? fSDigitsManager?
225     delete fSDigitsManager;
226     fSDigitsManager     = 0;
227   }
228
229   if (fSDigitsManagerList) {
230     fSDigitsManagerList->Delete();
231     delete fSDigitsManagerList;
232     fSDigitsManagerList = 0;
233   }
234
235   if (fMasks) {
236     delete [] fMasks;
237     fMasks       = 0;
238   }
239
240   if (fGeo) {
241     delete fGeo;
242     fGeo = 0;
243   }
244
245 }
246
247 //_____________________________________________________________________________
248 AliTRDdigitizer &AliTRDdigitizer::operator=(const AliTRDdigitizer &d)
249 {
250   //
251   // Assignment operator
252   //
253
254   if (this != &d) {
255     ((AliTRDdigitizer &) d).Copy(*this);
256   }
257
258   return *this;
259
260 }
261
262 //_____________________________________________________________________________
263 void AliTRDdigitizer::Copy(TObject &d) const
264 {
265   //
266   // Copy function
267   //
268
269   ((AliTRDdigitizer &) d).fRunLoader          = 0;
270   ((AliTRDdigitizer &) d).fDigitsManager      = 0;
271   ((AliTRDdigitizer &) d).fSDigitsManager     = 0;
272   ((AliTRDdigitizer &) d).fSDigitsManagerList = 0;
273   ((AliTRDdigitizer &) d).fTRD                = 0;
274   ((AliTRDdigitizer &) d).fGeo                = 0;
275   ((AliTRDdigitizer &) d).fEvent              = 0;
276   ((AliTRDdigitizer &) d).fMasks              = 0;
277   ((AliTRDdigitizer &) d).fCompress           = fCompress;
278   ((AliTRDdigitizer &) d).fSDigits            = fSDigits;
279   ((AliTRDdigitizer &) d).fMergeSignalOnly    = fMergeSignalOnly;
280
281 }
282
283 //_____________________________________________________________________________
284 void AliTRDdigitizer::Exec(const Option_t * const option)
285 {
286   //
287   // Executes the merging
288   //
289
290   Int_t iInput;
291
292   AliTRDdigitsManager *sdigitsManager;
293
294   TString optionString = option;
295   if (optionString.Contains("deb")) {
296     AliLog::SetClassDebugLevel("AliTRDdigitizer",1);
297     AliInfo("Called with debug option");
298   }
299
300   // The AliRoot file is already connected by the manager
301   AliRunLoader *inrl = 0x0;
302   
303   if (gAlice) {
304     AliDebug(1,"AliRun object found on file.");
305   }
306   else {
307     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
308     inrl->LoadgAlice();
309     gAlice = inrl->GetAliRun();
310     if (!gAlice) {
311       AliError("Could not find AliRun object.")
312       return;
313     }
314   }
315                                                                            
316   Int_t nInput = fManager->GetNinputs();
317   fMasks       = new Int_t[nInput];
318   for (iInput = 0; iInput < nInput; iInput++) {
319     fMasks[iInput] = fManager->GetMask(iInput);
320   }
321
322   //
323   // Initialization
324   //
325
326   AliRunLoader *orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
327
328   if (InitDetector()) {
329
330     AliLoader *ogime = orl->GetLoader("TRDLoader");
331
332     TTree *tree = 0;
333     if (fSDigits) { 
334       // If we produce SDigits
335       tree = ogime->TreeS();
336       if (!tree) {
337         ogime->MakeTree("S");
338         tree = ogime->TreeS();
339       }
340     }
341     else {
342       // If we produce Digits
343       tree = ogime->TreeD();
344       if (!tree) {
345         ogime->MakeTree("D");
346         tree = ogime->TreeD();
347       }
348     }
349
350     MakeBranch(tree);
351
352   }
353  
354   for (iInput = 0; iInput < nInput; iInput++) {
355
356     AliDebug(1,Form("Add input stream %d",iInput));
357
358     // Check if the input tree exists
359     inrl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
360     AliLoader *gime = inrl->GetLoader("TRDLoader");
361
362     TTree *treees = gime->TreeS();
363     if (treees == 0x0) {
364       if (gime->LoadSDigits()) {
365         AliError(Form("Error Occured while loading S. Digits for input %d.",iInput));
366         return;
367       }
368       treees = gime->TreeS();
369     }
370     
371     if (treees == 0x0) {
372       AliError(Form("Input stream %d does not exist",iInput));
373       return;
374     } 
375
376     // Read the s-digits via digits manager
377     sdigitsManager = new AliTRDdigitsManager();
378     sdigitsManager->SetSDigits(kTRUE);
379     
380     AliRunLoader *rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
381     AliLoader *gimme = rl->GetLoader("TRDLoader");
382     if (!gimme->TreeS()) 
383       {
384         gimme->LoadSDigits();
385       }
386
387     sdigitsManager->ReadDigits(gimme->TreeS());
388    
389     // Add the s-digits to the input list 
390     AddSDigitsManager(sdigitsManager);
391
392   }
393
394   // Convert the s-digits to normal digits
395   AliDebug(1,"Do the conversion");
396   SDigits2Digits();
397
398   // Store the digits
399   AliDebug(1,"Write the digits");
400   fDigitsManager->WriteDigits();
401
402   // Write parameters
403   orl->CdGAFile();
404
405   // Clean up
406   DeleteSDigitsManager();
407
408   AliDebug(1,"Done");
409
410 }
411
412 //_____________________________________________________________________________
413 Bool_t AliTRDdigitizer::Open(const Char_t *file, Int_t nEvent)
414 {
415   //
416   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
417   //
418   // Connect the AliRoot file containing Geometry, Kine, and Hits
419   //  
420
421   TString evfoldname = AliConfig::GetDefaultEventFolderName();
422
423   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
424   if (!fRunLoader) {
425     fRunLoader = AliRunLoader::Open(file,evfoldname,"UPDATE");
426   }  
427   if (!fRunLoader) {
428     AliError(Form("Can not open session for file %s.",file));
429     return kFALSE;
430   }
431    
432   if (!fRunLoader->GetAliRun()) {
433     fRunLoader->LoadgAlice();
434   }
435   gAlice = fRunLoader->GetAliRun();
436   
437   if (gAlice) {
438     AliDebug(1,"AliRun object found on file.");
439   }
440   else {
441     AliError("Could not find AliRun object.");
442     return kFALSE;
443   }
444
445   fEvent = nEvent;
446
447   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
448   if (!loader) {
449     AliError("Can not get TRD loader from Run Loader");
450     return kFALSE;
451   }
452   
453   if (InitDetector()) {
454     TTree *tree = 0;
455     if (fSDigits) { 
456       // If we produce SDigits
457       tree = loader->TreeS();
458       if (!tree) {
459         loader->MakeTree("S");
460         tree = loader->TreeS();
461       }
462     }
463     else {
464       // If we produce Digits
465       tree = loader->TreeD();
466       if (!tree) {
467         loader->MakeTree("D");
468         tree = loader->TreeD();
469       }
470     }
471     return MakeBranch(tree);
472   }
473   else {
474     return kFALSE;
475   }
476
477 }
478
479 //_____________________________________________________________________________
480 Bool_t AliTRDdigitizer::Open(AliRunLoader * const runLoader, Int_t nEvent)
481 {
482   //
483   // Opens a ROOT-file with TRD-hits and reads in the hit-tree
484   //
485   // Connect the AliRoot file containing Geometry, Kine, and Hits
486   //  
487
488   fRunLoader = runLoader;
489   if (!fRunLoader) {
490     AliError("RunLoader does not exist");
491     return kFALSE;
492   }
493    
494   if (!fRunLoader->GetAliRun()) {
495     fRunLoader->LoadgAlice();
496   }
497   gAlice = fRunLoader->GetAliRun();
498   
499   if (gAlice) {
500     AliDebug(1,"AliRun object found on file.");
501   }
502   else {
503     AliError("Could not find AliRun object.");
504     return kFALSE;
505   }
506
507   fEvent = nEvent;
508
509   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
510   if (!loader) {
511     AliError("Can not get TRD loader from Run Loader");
512     return kFALSE;
513   }
514   
515   if (InitDetector()) {
516     TTree *tree = 0;
517     if (fSDigits) { 
518       // If we produce SDigits
519       tree = loader->TreeS();
520       if (!tree) {
521         loader->MakeTree("S");
522         tree = loader->TreeS();
523       }
524     }
525     else {
526       // If we produce Digits
527       tree = loader->TreeD();
528       if (!tree) {
529         loader->MakeTree("D");
530         tree = loader->TreeD();
531       }
532     }
533     return MakeBranch(tree);
534   }
535   else {
536     return kFALSE;
537   }
538
539 }
540
541 //_____________________________________________________________________________
542 Bool_t AliTRDdigitizer::InitDetector()
543 {
544   //
545   // Sets the pointer to the TRD detector and the geometry
546   //
547
548   // Get the pointer to the detector class and check for version 1
549   fTRD = (AliTRD *) gAlice->GetDetector("TRD");
550   if (!fTRD) {
551     AliFatal("No TRD module found");
552     exit(1);
553   }
554   if (fTRD->IsVersion() != 1) {
555     AliFatal("TRD must be version 1 (slow simulator)");
556     exit(1);
557   }
558
559   // Get the geometry
560   fGeo = new AliTRDgeometry();
561
562   // Create a digits manager
563   if (fDigitsManager) {
564     delete fDigitsManager;
565   }
566   fDigitsManager = new AliTRDdigitsManager();
567   fDigitsManager->SetSDigits(fSDigits);
568   fDigitsManager->CreateArrays();
569   fDigitsManager->SetEvent(fEvent);
570
571   // The list for the input s-digits manager to be merged
572   if (fSDigitsManagerList) {
573     fSDigitsManagerList->Delete();
574   } 
575   else {
576     fSDigitsManagerList = new TList();
577   }
578
579   return kTRUE;
580
581 }
582 //_____________________________________________________________________________
583 Bool_t AliTRDdigitizer::MakeBranch(TTree *tree) const
584 {
585   // 
586   // Create the branches for the digits array
587   //
588
589   return fDigitsManager->MakeBranch(tree);
590
591 }
592
593 //_____________________________________________________________________________
594 void AliTRDdigitizer::AddSDigitsManager(AliTRDdigitsManager *man)
595 {
596   //
597   // Add a digits manager for s-digits to the input list.
598   //
599
600   fSDigitsManagerList->Add(man);
601
602 }
603
604 //_____________________________________________________________________________
605 void AliTRDdigitizer::DeleteSDigitsManager()
606 {
607   //
608   // Removes digits manager from the input list.
609   //
610
611   fSDigitsManagerList->Delete();
612
613 }
614
615 //_____________________________________________________________________________
616 Bool_t AliTRDdigitizer::MakeDigits()
617 {
618   //
619   // Creates digits.
620   //
621
622   AliDebug(1,"Start creating digits");
623
624   if (!fGeo) {
625     AliError("No geometry defined");
626     return kFALSE;
627   }
628
629   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
630   if (!calibration) {
631     AliFatal("Could not get calibration object");
632     return kFALSE;
633   }
634
635   const Int_t kNdet  = AliTRDgeometry::Ndet();
636
637   Float_t **hits = new Float_t*[kNdet];
638   Int_t    *nhit = new Int_t[kNdet];
639
640   AliTRDarraySignal *signals = 0x0;
641  
642   // Sort all hits according to detector number
643   if (!SortHits(hits,nhit)) {
644     AliError("Sorting hits failed");
645     return kFALSE;
646   }
647
648   // Loop through all detectors
649   for (Int_t det = 0; det < kNdet; det++) {
650
651     // Detectors that are switched off, not installed, etc.
652     if (( calibration->IsChamberInstalled(det)) &&
653         (!calibration->IsChamberMasked(det))    &&
654         ( fGeo->ChamberInGeometry(det))         &&
655         (nhit[det] > 0)) {
656
657       signals = new AliTRDarraySignal();
658           
659       // Convert the hits of the current detector to detector signals
660       if (!ConvertHits(det,hits[det],nhit[det],signals)) {
661         AliError(Form("Conversion of hits failed for detector=%d",det));
662         return kFALSE;
663       }
664       // Convert the detector signals to digits or s-digits
665       if (!ConvertSignals(det,signals)) {
666         AliError(Form("Conversion of signals failed for detector=%d",det));
667         return kFALSE;
668       }
669
670       // Delete the signals array
671       delete signals;
672       signals = 0x0;
673
674     } // if: detector status
675
676     delete [] hits[det];
677
678   } // for: detector
679
680   delete [] hits;
681   delete [] nhit;
682
683   return kTRUE;
684
685 }
686
687 //_____________________________________________________________________________
688 Bool_t AliTRDdigitizer::SortHits(Float_t **hits, Int_t *nhit)
689 {
690   //
691   // Read all the hits and sorts them according to detector number
692   // in the output array <hits>.
693   //
694
695   AliDebug(1,"Start sorting hits");
696
697   const Int_t kNdet = AliTRDgeometry::Ndet();
698   // Size of the hit vector
699   const Int_t kNhit = 6;
700
701   Float_t *xyz      = 0;
702   Int_t    nhitTrk  = 0;
703
704   Int_t   *lhit     = new Int_t[kNdet];
705
706   for (Int_t det = 0; det < kNdet; det++) {
707     lhit[det] = 0;
708     nhit[det] = 0;
709     hits[det] = 0;
710   }
711
712   AliLoader *gimme = fRunLoader->GetLoader("TRDLoader");
713   if (!gimme->TreeH()) {
714     gimme->LoadHits();
715   }
716   TTree *hitTree = gimme->TreeH();
717   if (hitTree == 0x0) {
718     AliError("Can not get TreeH");
719     return kFALSE;
720   }
721   fTRD->SetTreeAddress();
722
723   // Get the number of entries in the hit tree
724   // (Number of primary particles creating a hit somewhere)
725   Int_t nTrk = (Int_t) hitTree->GetEntries();
726   AliDebug(1,Form("Found %d tracks",nTrk));
727
728   // Loop through all the tracks in the tree
729   for (Int_t iTrk = 0; iTrk < nTrk; iTrk++) {
730
731     gAlice->GetMCApp()->ResetHits();
732     hitTree->GetEvent(iTrk);
733
734     if (!fTRD->Hits()) {
735       AliError(Form("No hits array for track = %d",iTrk));
736       continue;
737     }
738
739     // Number of hits for this track
740     nhitTrk = fTRD->Hits()->GetEntriesFast();
741
742     Int_t hitCnt = 0;
743     // Loop through the TRD hits
744     AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1);
745     while (hit) {
746
747       hitCnt++;
748
749       // Don't analyze test hits
750       if (((Int_t) hit->GetCharge()) != 0) {
751
752         Int_t   trk  = hit->Track();
753         Int_t   det  = hit->GetDetector();
754         Int_t   q    = hit->GetCharge();
755         Float_t x    = hit->X();
756         Float_t y    = hit->Y();
757         Float_t z    = hit->Z();
758         Float_t time = hit->GetTime();
759
760         if (nhit[det] == lhit[det]) {
761           // Inititialization of new detector
762           xyz        = new Float_t[kNhit*(nhitTrk+lhit[det])];
763           if (hits[det]) {
764             memcpy(xyz,hits[det],sizeof(Float_t)*kNhit*lhit[det]);
765             delete [] hits[det];
766           }
767           lhit[det] += nhitTrk;
768           hits[det]  = xyz;
769         }
770         else {
771           xyz        = hits[det];
772         }
773         xyz[nhit[det]*kNhit+0] = x;
774         xyz[nhit[det]*kNhit+1] = y;
775         xyz[nhit[det]*kNhit+2] = z;
776         xyz[nhit[det]*kNhit+3] = q;
777         xyz[nhit[det]*kNhit+4] = trk;
778         xyz[nhit[det]*kNhit+5] = time;
779         nhit[det]++;
780
781       } // if: charge != 0
782
783       hit = (AliTRDhit *) fTRD->NextHit();   
784
785     } // for: hits of one track
786
787   } // for: tracks
788
789   delete [] lhit;
790
791   return kTRUE;
792
793 }
794
795 //_____________________________________________________________________________
796 Bool_t AliTRDdigitizer::ConvertHits(Int_t det
797                                   , const Float_t * const hits
798                                   , Int_t nhit
799                                   , AliTRDarraySignal *signals)
800 {
801   //
802   // Converts the detectorwise sorted hits to detector signals
803   //
804
805   AliDebug(1,Form("Start converting hits for detector=%d (nhits=%d)",det,nhit));
806
807   // Number of pads included in the pad response
808   const Int_t kNpad      = 3;
809   // Number of track dictionary arrays
810   const Int_t kNdict     = AliTRDdigitsManager::kNDict;
811   // Size of the hit vector
812   const Int_t kNhit      = 6;
813
814   // Width of the amplification region
815   const Float_t kAmWidth = AliTRDgeometry::AmThick();
816   // Width of the drift region
817   const Float_t kDrWidth = AliTRDgeometry::DrThick();
818   // Drift + amplification region 
819   const Float_t kDrMin   =          - 0.5 * kAmWidth;
820   const Float_t kDrMax   = kDrWidth + 0.5 * kAmWidth;
821   
822   Int_t    iPad          = 0;
823   Int_t    dict          = 0;
824   Int_t    timeBinTRFend = 1;
825
826   Double_t pos[3];
827   Double_t loc[3];
828   Double_t padSignal[kNpad];
829   Double_t signalOld[kNpad];
830
831   AliTRDarrayDictionary *dictionary[kNdict];
832
833   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
834   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
835   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
836
837   if (!commonParam) {
838     AliFatal("Could not get common parameterss");
839     return kFALSE;
840   }
841   if (!simParam) {
842     AliFatal("Could not get simulation parameters");
843     return kFALSE;
844   }  
845   if (!calibration) {
846     AliFatal("Could not get calibration object");  
847     return kFALSE;
848   }
849
850   // Get the detector wise calibration objects
851   AliTRDCalROC       *calVdriftROC      = 0;
852   Float_t             calVdriftDetValue = 0.0;
853   const AliTRDCalDet *calVdriftDet      = calibration->GetVdriftDet();  
854   AliTRDCalROC       *calT0ROC          = 0;
855   Float_t             calT0DetValue     = 0.0;
856   const AliTRDCalDet *calT0Det          = calibration->GetT0Det();  
857
858   if (simParam->TRFOn()) {
859     timeBinTRFend = ((Int_t) (simParam->GetTRFhi() 
860                   * commonParam->GetSamplingFrequency())) - 1;
861   }
862
863   Int_t   nTimeTotal   = calibration->GetNumberOfTimeBins();
864   Float_t samplingRate = commonParam->GetSamplingFrequency();
865   Float_t elAttachProp = simParam->GetElAttachProp() / 100.0; 
866
867   AliTRDpadPlane *padPlane = fGeo->GetPadPlane(det);
868   Int_t   layer   = fGeo->GetLayer(det);         //update
869   Float_t row0    = padPlane->GetRow0ROC();
870   Int_t   nRowMax = padPlane->GetNrows();
871   Int_t   nColMax = padPlane->GetNcols();
872
873   // Create a new array for the signals
874   signals->Allocate(nRowMax,nColMax,nTimeTotal);
875
876   // Create a new array for the dictionary
877   for (dict = 0; dict < kNdict; dict++) {       
878     dictionary[dict] = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
879     dictionary[dict]->Allocate(nRowMax,nColMax,nTimeTotal);
880   }      
881
882   // Loop through the hits in this detector
883   for (Int_t hit = 0; hit < nhit; hit++) {
884
885     pos[0]          = hits[hit*kNhit+0];
886     pos[1]          = hits[hit*kNhit+1];
887     pos[2]          = hits[hit*kNhit+2];
888     Float_t q       = hits[hit*kNhit+3];
889     Float_t hittime = hits[hit*kNhit+5];
890     Int_t   track   = ((Int_t) hits[hit*kNhit+4]);
891
892     Int_t   inDrift = 1;
893
894     // Find the current volume with the geo manager
895     gGeoManager->SetCurrentPoint(pos);
896     gGeoManager->FindNode();      
897     if (strstr(gGeoManager->GetPath(),"/UK")) {
898       inDrift = 0;
899     }
900
901     // Get the calibration objects
902     calVdriftROC      = calibration->GetVdriftROC(det);
903     calVdriftDetValue = calVdriftDet->GetValue(det);
904     calT0ROC          = calibration->GetT0ROC(det);
905     calT0DetValue     = calT0Det->GetValue(det);
906
907     // Go to the local coordinate system:
908     // loc[0] - col  direction in amplification or driftvolume
909     // loc[1] - row  direction in amplification or driftvolume
910     // loc[2] - time direction in amplification or driftvolume
911     gGeoManager->MasterToLocal(pos,loc);
912     if (inDrift) {
913       // Relative to middle of amplification region
914       loc[2] = loc[2] - kDrWidth/2.0 - kAmWidth/2.0;
915     } 
916
917     // The driftlength [cm] (w/o diffusion yet !).
918     // It is negative if the hit is between pad plane and anode wires.
919     Double_t driftlength = -1.0 * loc[2];
920
921     // Stupid patch to take care of TR photons that are absorbed
922     // outside the chamber volume. A real fix would actually need
923     // a more clever implementation of the TR hit generation
924     if (q < 0.0) {
925       if ((loc[1] < padPlane->GetRowEndROC()) ||
926           (loc[1] > padPlane->GetRow0ROC())) {
927         continue;
928       }
929       if ((driftlength < kDrMin) ||
930           (driftlength > kDrMax)) {
931         continue;
932       }
933     }
934
935     // Get row and col of unsmeared electron to retrieve drift velocity
936     // The pad row (z-direction)
937     Int_t    rowE         = padPlane->GetPadRowNumberROC(loc[1]);
938     if (rowE < 0) {
939       continue;
940     }
941     Double_t rowOffset    = padPlane->GetPadRowOffsetROC(rowE,loc[1]);
942     // The pad column (rphi-direction)
943     Double_t offsetTilt   = padPlane->GetTiltOffset(rowOffset);
944     Int_t    colE         = padPlane->GetPadColNumber(loc[0]+offsetTilt);
945     if (colE < 0) {
946       continue;   
947     }
948     Double_t colOffset    = 0.0;
949
950     // Normalized drift length
951     Float_t  driftvelocity  = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
952     Double_t absdriftlength = TMath::Abs(driftlength);
953     if (commonParam->ExBOn()) {
954       absdriftlength /= TMath::Sqrt(GetLorentzFactor(driftvelocity));
955     }
956
957     // Loop over all electrons of this hit
958     // TR photons produce hits with negative charge
959     Int_t nEl = ((Int_t) TMath::Abs(q));
960     for (Int_t iEl = 0; iEl < nEl; iEl++) {
961
962       // Now the real local coordinate system of the ROC
963       // column direction: locC
964       // row direction:    locR 
965       // time direction:   locT
966       // locR and locC are identical to the coordinates of the corresponding
967       // volumina of the drift or amplification region.
968       // locT is defined relative to the wire plane (i.e. middle of amplification
969       // region), meaning locT = 0, and is negative for hits coming from the
970       // drift region. 
971       Double_t locC = loc[0];
972       Double_t locR = loc[1];
973       Double_t locT = loc[2];
974
975       // Electron attachment
976       if (simParam->ElAttachOn()) {
977         if (gRandom->Rndm() < (absdriftlength * elAttachProp)) {
978           continue;
979         }
980       }
981           
982       // Apply the diffusion smearing
983       if (simParam->DiffusionOn()) {
984         if (!(Diffusion(driftvelocity,absdriftlength,locR,locC,locT))) {
985           continue;
986         }
987       }
988
989       // Apply E x B effects (depends on drift direction)
990       if (commonParam->ExBOn()) { 
991         if (!(ExB(driftvelocity,driftlength,locC))) {
992           continue;
993         }
994       }
995
996       // The electron position after diffusion and ExB in pad coordinates.
997       // The pad row (z-direction)
998       rowE       = padPlane->GetPadRowNumberROC(locR);
999       if (rowE < 0) continue;
1000       rowOffset  = padPlane->GetPadRowOffsetROC(rowE,locR);
1001
1002       // The pad column (rphi-direction)
1003       offsetTilt = padPlane->GetTiltOffset(rowOffset);
1004       colE       = padPlane->GetPadColNumber(locC+offsetTilt);
1005       if (colE < 0) continue;         
1006       colOffset  = padPlane->GetPadColOffset(colE,locC+offsetTilt);
1007           
1008       // Also re-retrieve drift velocity because col and row may have changed
1009       driftvelocity = calVdriftDetValue * calVdriftROC->GetValue(colE,rowE);
1010       Float_t t0    = calT0DetValue     + calT0ROC->GetValue(colE,rowE);
1011
1012       // Convert the position to drift time [mus], using either constant drift velocity or
1013       // time structure of drift cells (non-isochronity, GARFIELD calculation).
1014       // Also add absolute time of hits to take pile-up events into account properly
1015       Double_t drifttime;
1016       if (simParam->TimeStructOn()) {
1017         // Get z-position with respect to anode wire
1018         Double_t zz  =  row0 - locR + padPlane->GetAnodeWireOffset();
1019         zz -= ((Int_t)(2 * zz)) / 2.0;
1020         if (zz > 0.25) {
1021           zz  = 0.5 - zz;
1022         }
1023         // Use drift time map (GARFIELD)
1024         drifttime = commonParam->TimeStruct(driftvelocity,0.5*kAmWidth-1.0*locT,zz)
1025                   + hittime;
1026       } 
1027       else {
1028         // Use constant drift velocity
1029         drifttime = TMath::Abs(locT) / driftvelocity
1030                   + hittime;
1031       }
1032
1033       // Apply the gas gain including fluctuations
1034       Double_t ggRndm = 0.0;
1035       do {
1036         ggRndm = gRandom->Rndm();
1037       } while (ggRndm <= 0);
1038       Double_t signal = -(simParam->GetGasGain()) * TMath::Log(ggRndm);
1039
1040       // Apply the pad response 
1041       if (simParam->PRFOn()) {
1042         // The distance of the electron to the center of the pad 
1043         // in units of pad width
1044         Double_t dist = (colOffset - 0.5*padPlane->GetColSize(colE))
1045                       / padPlane->GetColSize(colE);
1046         // This is a fixed parametrization, i.e. not dependent on
1047         // calibration values !
1048         if (!(calibration->PadResponse(signal,dist,layer,padSignal))) continue;
1049       }
1050       else {
1051         padSignal[0] = 0.0;
1052         padSignal[1] = signal;
1053         padSignal[2] = 0.0;
1054       }
1055
1056       // The time bin (always positive), with t0 distortion
1057       Double_t timeBinIdeal = drifttime * samplingRate + t0;
1058       // Protection 
1059       if (TMath::Abs(timeBinIdeal) > 2*nTimeTotal) {
1060         timeBinIdeal = 2 * nTimeTotal;
1061       }
1062       Int_t    timeBinTruncated = ((Int_t) timeBinIdeal);
1063       // The distance of the position to the middle of the timebin
1064       Double_t timeOffset       = ((Float_t) timeBinTruncated 
1065                                 + 0.5 - timeBinIdeal) / samplingRate;
1066           
1067       // Sample the time response inside the drift region
1068       // + additional time bins before and after.
1069       // The sampling is done always in the middle of the time bin
1070       for (Int_t iTimeBin = TMath::Max(timeBinTruncated,0)
1071           ;iTimeBin < TMath::Min(timeBinTruncated+timeBinTRFend,nTimeTotal)
1072           ;iTimeBin++) {
1073
1074         // Apply the time response
1075         Double_t timeResponse = 1.0;
1076         Double_t crossTalk    = 0.0;
1077         Double_t time         = (iTimeBin - timeBinTruncated) / samplingRate + timeOffset;
1078
1079         if (simParam->TRFOn()) {
1080           timeResponse = simParam->TimeResponse(time);
1081         }
1082         if (simParam->CTOn()) {
1083           crossTalk    = simParam->CrossTalk(time);
1084         }
1085
1086         signalOld[0] = 0.0;
1087         signalOld[1] = 0.0;
1088         signalOld[2] = 0.0;
1089
1090         for (iPad = 0; iPad < kNpad; iPad++) {
1091
1092           Int_t colPos = colE + iPad - 1;
1093           if (colPos <        0) continue;
1094           if (colPos >= nColMax) break;
1095
1096           // Add the signals
1097           signalOld[iPad] = signals->GetData(rowE,colPos,iTimeBin);
1098
1099           if (colPos != colE) {
1100             // Cross talk added to non-central pads
1101             signalOld[iPad] += padSignal[iPad] 
1102                              * (timeResponse + crossTalk);
1103           } 
1104           else {
1105             // W/o cross talk at central pad
1106             signalOld[iPad] += padSignal[iPad] 
1107                              * timeResponse;
1108           }
1109
1110           signals->SetData(rowE,colPos,iTimeBin,signalOld[iPad]);
1111
1112           // Store the track index in the dictionary
1113           // Note: We store index+1 in order to allow the array to be compressed
1114           // Note2: Taking out the +1 in track
1115           if (signalOld[iPad] > 0.0) { 
1116             for (dict = 0; dict < kNdict; dict++) {
1117               Int_t oldTrack = dictionary[dict]->GetData(rowE,colPos,iTimeBin); 
1118               if (oldTrack == track) break;
1119               if (oldTrack ==  -1 ) {
1120                 dictionary[dict]->SetData(rowE,colPos,iTimeBin,track);
1121                 break;
1122               }
1123             }
1124           }
1125
1126         } // Loop: pads
1127
1128       } // Loop: time bins
1129
1130     } // Loop: electrons of a single hit
1131
1132   } // Loop: hits
1133
1134   AliDebug(2,Form("Finished analyzing %d hits",nhit));
1135
1136   return kTRUE;
1137
1138 }
1139
1140 //_____________________________________________________________________________
1141 Bool_t AliTRDdigitizer::ConvertSignals(Int_t det, AliTRDarraySignal *signals)
1142 {
1143   //
1144   // Convert signals to digits
1145   //
1146
1147   AliDebug(1,Form("Start converting the signals for detector %d",det));
1148
1149   if (fSDigits) {
1150     // Convert the signal array to s-digits
1151     if (!Signal2SDigits(det,signals)) {
1152       return kFALSE;
1153     }
1154   }
1155   else {
1156     // Convert the signal array to digits
1157     if (!Signal2ADC(det,signals)) {
1158       return kFALSE;
1159     }
1160     // Run digital processing for digits
1161     RunDigitalProcessing(det);
1162   }
1163
1164   // Compress the arrays
1165   CompressOutputArrays(det);   
1166
1167   return kTRUE;
1168
1169 }
1170
1171 //_____________________________________________________________________________
1172 Bool_t AliTRDdigitizer::Signal2ADC(Int_t det, AliTRDarraySignal *signals)
1173 {
1174   //
1175   // Converts the sampled electron signals to ADC values for a given chamber
1176   //
1177
1178   AliDebug(1,Form("Start converting signals to ADC values for detector=%d",det));
1179
1180   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1181   if (!calibration) {
1182     AliFatal("Could not get calibration object");
1183     return kFALSE;
1184   }
1185
1186   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1187   if (!simParam) {
1188     AliFatal("Could not get simulation parameters");
1189     return kFALSE;
1190   }  
1191
1192   // Converts number of electrons to fC
1193   const Double_t kEl2fC = 1.602e-19 * 1.0e15; 
1194
1195   // Coupling factor
1196   Double_t coupling     = simParam->GetPadCoupling() 
1197                         * simParam->GetTimeCoupling();
1198   // Electronics conversion factor
1199   Double_t convert      = kEl2fC 
1200                         * simParam->GetChipGain();
1201   // ADC conversion factor
1202   Double_t adcConvert   = simParam->GetADCoutRange()
1203                         / simParam->GetADCinRange();
1204   // The electronics baseline in mV
1205   Double_t baseline     = simParam->GetADCbaseline() 
1206                         / adcConvert;
1207   // The electronics baseline in electrons
1208   Double_t baselineEl   = baseline
1209                         / convert;
1210
1211   Int_t row  = 0;
1212   Int_t col  = 0;
1213   Int_t time = 0;
1214
1215   Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1216   Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1217   Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
1218
1219   // The gainfactor calibration objects
1220   const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();  
1221   AliTRDCalROC       *calGainFactorROC      = 0;
1222   Float_t             calGainFactorDetValue = 0.0;
1223
1224   AliTRDarrayADC     *digits = 0x0;
1225
1226   if (!signals) {
1227     AliError(Form("Signals array for detector %d does not exist\n",det));
1228     return kFALSE;
1229   }
1230   if (signals->HasData()) {
1231     // Expand the container if neccessary
1232     signals->Expand();   
1233   }
1234   else {
1235     // Create missing containers
1236     signals->Allocate(nRowMax,nColMax,nTimeTotal);      
1237   }
1238
1239   // Get the container for the digits of this detector
1240   if (fDigitsManager->HasSDigits()) {
1241     AliError("Digits manager has s-digits");
1242     return kFALSE;
1243   }
1244
1245   digits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det);
1246   // Allocate memory space for the digits buffer
1247   if (!digits->HasData()) {
1248     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1249   }
1250
1251   // Get the calibration objects
1252   calGainFactorROC      = calibration->GetGainFactorROC(det);
1253   calGainFactorDetValue = calGainFactorDet->GetValue(det);
1254
1255   // Create the digits for this chamber
1256   for (row  = 0; row  <  nRowMax; row++ ) {
1257     for (col  = 0; col  <  nColMax; col++ ) {
1258
1259       // Check whether pad is masked
1260       // Bridged pads are not considered yet!!!
1261       if (calibration->IsPadMasked(det,col,row)) {
1262         continue;
1263       }
1264
1265       // The gain factors
1266       Float_t padgain = calGainFactorDetValue 
1267                       * calGainFactorROC->GetValue(col,row);
1268       if (padgain <= 0) {
1269         AliError(Form("Not a valid gain %f, %d %d %d",padgain,det,col,row));
1270       }
1271
1272       for (time = 0; time < nTimeTotal; time++) {
1273
1274         // Get the signal amplitude
1275         Float_t signalAmp = signals->GetData(row,col,time);
1276         // Pad and time coupling
1277         signalAmp *= coupling;
1278         // Gain factors
1279         signalAmp *= padgain;
1280
1281         // Add the noise, starting from minus ADC baseline in electrons
1282         signalAmp  = TMath::Max((Double_t) gRandom->Gaus(signalAmp,simParam->GetNoise())
1283                                ,-baselineEl);
1284
1285         // Convert to mV
1286         signalAmp *= convert;
1287         // Add ADC baseline in mV
1288         signalAmp += baseline;
1289
1290         // Convert to ADC counts. Set the overflow-bit fADCoutRange if the
1291         // signal is larger than fADCinRange
1292         Short_t adc  = 0;
1293         if (signalAmp >= simParam->GetADCinRange()) {
1294           adc = ((Short_t) simParam->GetADCoutRange());
1295         }
1296         else {
1297           adc = TMath::Nint(signalAmp * adcConvert);
1298         }
1299
1300         // Saving all digits
1301         digits->SetData(row,col,time,adc);
1302
1303       } // for: time
1304
1305     } // for: col
1306   } // for: row
1307
1308   return kTRUE;
1309
1310 }
1311
1312 //_____________________________________________________________________________
1313 Bool_t AliTRDdigitizer::Signal2SDigits(Int_t det, AliTRDarraySignal *signals)
1314 {
1315   //
1316   // Converts the sampled electron signals to s-digits
1317   //
1318
1319   AliDebug(1,Form("Start converting signals to s-digits for detector=%d",det));
1320
1321   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1322   if (!calibration) {
1323     AliFatal("Could not get calibration object");
1324     return kFALSE;
1325   }
1326
1327   Int_t row  = 0;
1328   Int_t col  = 0;
1329   Int_t time = 0;
1330
1331   Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1332   Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1333   Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
1334
1335   // Get the container for the digits of this detector
1336
1337   if (!fDigitsManager->HasSDigits()) {
1338     AliError("Digits manager has no s-digits");
1339     return kFALSE;
1340   }
1341
1342   AliTRDarraySignal *digits = (AliTRDarraySignal *) fDigitsManager->GetSDigits(det);
1343   // Allocate memory space for the digits buffer
1344   if (!digits->HasData()) {
1345     digits->Allocate(nRowMax,nColMax,nTimeTotal);
1346   }
1347
1348   // Create the sdigits for this chamber
1349   for (row  = 0; row  <  nRowMax; row++ ) {
1350     for (col  = 0; col  <  nColMax; col++ ) {
1351       for (time = 0; time < nTimeTotal; time++) {         
1352         digits->SetData(row,col,time,signals->GetData(row,col,time));
1353       } // for: time
1354     } // for: col
1355   } // for: row
1356   
1357   return kTRUE;
1358
1359 }
1360
1361 //_____________________________________________________________________________
1362 Bool_t AliTRDdigitizer::Digits2SDigits(AliTRDdigitsManager * const manDig
1363                                      , AliTRDdigitsManager * const manSDig)
1364 {
1365   //
1366   // Converts digits into s-digits. Needed for embedding into real data.
1367   //
1368
1369   AliDebug(1,"Start converting digits to s-digits");
1370
1371   if (!fGeo) {
1372     fGeo = new AliTRDgeometry();
1373   }
1374
1375   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1376   if (!calibration) {
1377     AliFatal("Could not get calibration object");
1378     return kFALSE;
1379   }
1380
1381   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1382   if (!simParam) {
1383     AliFatal("Could not get simulation parameters");
1384     return kFALSE;
1385   }  
1386
1387   // Converts number of electrons to fC
1388   const Double_t kEl2fC = 1.602e-19 * 1.0e15; 
1389
1390   // Coupling factor
1391   Double_t coupling     = simParam->GetPadCoupling() 
1392                         * simParam->GetTimeCoupling();
1393   // Electronics conversion factor
1394   Double_t convert      = kEl2fC 
1395                         * simParam->GetChipGain();
1396   // ADC conversion factor
1397   Double_t adcConvert   = simParam->GetADCoutRange()
1398                         / simParam->GetADCinRange();
1399   // The electronics baseline in mV
1400   Double_t baseline     = simParam->GetADCbaseline() 
1401                         / adcConvert;
1402   // The electronics baseline in electrons
1403   //Double_t baselineEl   = baseline
1404   //                      / convert;
1405
1406   // The gainfactor calibration objects
1407   //const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();  
1408   //AliTRDCalROC       *calGainFactorROC      = 0;
1409   //Float_t             calGainFactorDetValue = 0.0;
1410
1411   Int_t row  = 0;
1412   Int_t col  = 0;
1413   Int_t time = 0;
1414
1415   for (Int_t det = 0; det < AliTRDgeometry::Ndet(); det++) {
1416
1417     Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1418     Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1419     Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
1420
1421     // Get the calibration objects
1422     //calGainFactorROC      = calibration->GetGainFactorROC(det);
1423     //calGainFactorDetValue = calGainFactorDet->GetValue(det);
1424
1425     // Get the digits
1426     AliTRDarrayADC *digits = (AliTRDarrayADC *) manDig->GetDigits(det);
1427
1428     if (!manSDig->HasSDigits()) {
1429       AliError("SDigits manager has no s-digits");
1430       return kFALSE;
1431     }
1432     // Get the s-digits
1433     AliTRDarraySignal     *sdigits = (AliTRDarraySignal *)     manSDig->GetSDigits(det);
1434     AliTRDarrayDictionary *tracks0 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,0);
1435     AliTRDarrayDictionary *tracks1 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,1);
1436     AliTRDarrayDictionary *tracks2 = (AliTRDarrayDictionary *) manSDig->GetDictionary(det,2);
1437     // Allocate memory space for the digits buffer
1438     sdigits->Allocate(nRowMax,nColMax,nTimeTotal);
1439     tracks0->Allocate(nRowMax,nColMax,nTimeTotal);
1440     tracks1->Allocate(nRowMax,nColMax,nTimeTotal);
1441     tracks2->Allocate(nRowMax,nColMax,nTimeTotal);
1442
1443     if (digits->HasData()) {
1444
1445       digits->Expand();
1446
1447       // Create the sdigits for this chamber
1448       for (row  = 0; row  <  nRowMax; row++ ) {
1449         for (col  = 0; col  <  nColMax; col++ ) {
1450
1451           // The gain factors
1452           //Float_t padgain = calGainFactorDetValue 
1453           //                * calGainFactorROC->GetValue(col,row);
1454
1455           for (time = 0; time < nTimeTotal; time++) {
1456
1457             Short_t  adcVal = digits->GetData(row,col,time);
1458             Double_t signal = (Double_t) adcVal;
1459             // ADC -> signal in mV
1460             signal /= adcConvert;
1461             // Subtract baseline in mV
1462             signal -= baseline;
1463             // Signal in mV -> signal in #electrons
1464             signal /= convert;
1465             // Gain factor
1466             //signal /= padgain; // Not needed for real data
1467             // Pad and time coupling
1468             signal /= coupling;
1469
1470             sdigits->SetData(row,col,time,signal);
1471             tracks0->SetData(row,col,time,0);
1472             tracks1->SetData(row,col,time,0);
1473             tracks2->SetData(row,col,time,0);
1474
1475           } // for: time
1476
1477         } // for: col
1478       } // for: row
1479   
1480     } // if: has data
1481
1482     sdigits->Compress(0);
1483     tracks0->Compress();
1484     tracks1->Compress();
1485     tracks2->Compress();
1486
1487     // No compress just remove
1488     manDig->RemoveDigits(det);
1489     manDig->RemoveDictionaries(det);      
1490
1491   } // for: det
1492
1493   return kTRUE;
1494
1495 }
1496
1497 //_____________________________________________________________________________
1498 Bool_t AliTRDdigitizer::SDigits2Digits()
1499 {
1500   //
1501   // Merges the input s-digits and converts them to normal digits
1502   //
1503
1504   if (!MergeSDigits()) {
1505     return kFALSE;
1506   }
1507
1508   return ConvertSDigits();
1509
1510 }
1511
1512 //_____________________________________________________________________________
1513 Bool_t AliTRDdigitizer::MergeSDigits()
1514 {
1515   //
1516   // Merges the input s-digits:
1517   //   - The amplitude of the different inputs are summed up.
1518   //   - Of the track IDs from the input dictionaries only one is
1519   //     kept for each input. This works for maximal 3 different merged inputs.
1520   //
1521
1522   // Number of track dictionary arrays
1523   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1524
1525   AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
1526   if (!simParam) {
1527     AliFatal("Could not get simulation parameters");
1528     return kFALSE;
1529   }
1530   
1531   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1532   if (!calibration) {
1533     AliFatal("Could not get calibration object");
1534     return kFALSE;
1535   }
1536   
1537   Int_t iDict = 0;
1538   Int_t jDict = 0;
1539
1540   AliTRDarraySignal     *digitsA;
1541   AliTRDarraySignal     *digitsB;
1542   AliTRDarrayDictionary *dictionaryA[kNDict];
1543   AliTRDarrayDictionary *dictionaryB[kNDict];
1544
1545   AliTRDdigitsManager   *mergeSDigitsManager = 0x0;
1546   // Get the first s-digits
1547   fSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->First();
1548   if (!fSDigitsManager) { 
1549     AliError("No SDigits manager");
1550     return kFALSE;
1551   }
1552
1553   // Loop through the other sets of s-digits
1554   mergeSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->After(fSDigitsManager);
1555
1556   if (mergeSDigitsManager) {
1557     AliDebug(1,Form("Merge %d input files.",fSDigitsManagerList->GetSize()));
1558   }
1559   else {
1560     AliDebug(1,"Only one input file.");
1561   }
1562   
1563   Int_t nTimeTotal = calibration->GetNumberOfTimeBins();  
1564   Int_t iMerge = 0;
1565
1566   while (mergeSDigitsManager) {
1567
1568     iMerge++;
1569       
1570     // Loop through the detectors
1571     for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) {
1572
1573       Int_t nRowMax = fGeo->GetPadPlane(iDet)->GetNrows();
1574       Int_t nColMax = fGeo->GetPadPlane(iDet)->GetNcols();
1575           
1576       // Loop through the pixels of one detector and add the signals
1577       digitsA = (AliTRDarraySignal *) fSDigitsManager->GetSDigits(iDet);    
1578       digitsB = (AliTRDarraySignal *) mergeSDigitsManager->GetSDigits(iDet); 
1579       digitsA->Expand();  
1580       if (!digitsA->HasData()) continue;
1581       digitsB->Expand();    
1582       if (!digitsB->HasData()) continue;
1583           
1584       for (iDict = 0; iDict < kNDict; iDict++) {
1585         dictionaryA[iDict] = (AliTRDarrayDictionary *) fSDigitsManager->GetDictionary(iDet,iDict);
1586         dictionaryB[iDict] = (AliTRDarrayDictionary *) mergeSDigitsManager->GetDictionary(iDet,iDict);
1587         dictionaryA[iDict]->Expand();  
1588         dictionaryB[iDict]->Expand();
1589       }
1590
1591       // Merge only detectors that contain a signal
1592       Bool_t doMerge = kTRUE;
1593       if (fMergeSignalOnly) {
1594         if (digitsA->GetOverThreshold(0) == 0) {                             
1595           doMerge = kFALSE;
1596         }
1597       }
1598           
1599       if (doMerge) {
1600               
1601         AliDebug(1,Form("Merge detector %d of input no.%d",iDet,iMerge+1));
1602               
1603         for (Int_t iRow  = 0; iRow  <  nRowMax;   iRow++ ) {
1604           for (Int_t iCol  = 0; iCol  <  nColMax;   iCol++ ) {
1605             for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
1606                 
1607               // Add the amplitudes of the summable digits 
1608               Float_t ampA = digitsA->GetData(iRow,iCol,iTime);
1609               Float_t ampB = digitsB->GetData(iRow,iCol,iTime);
1610               ampA += ampB;
1611               digitsA->SetData(iRow,iCol,iTime,ampA);
1612
1613               // Add the mask to the track id if defined.
1614               for (iDict = 0; iDict < kNDict; iDict++) {
1615                 Int_t trackB = dictionaryB[iDict]->GetData(iRow,iCol,iTime);
1616                 if ((fMasks) && (trackB > 0))  {
1617                   for (jDict = 0; jDict < kNDict; jDict++) { 
1618                     Int_t trackA = dictionaryA[iDict]->GetData(iRow,iCol,iTime); 
1619                     if (trackA == 0) {
1620                       trackA = trackB + fMasks[iMerge];
1621                       dictionaryA[iDict]->SetData(iRow,iCol,iTime,trackA);  
1622                     } // if:  track A == 0
1623                   } // for: jDict
1624                 } // if:  fMasks and trackB > 0
1625               } // for: iDict
1626
1627             } // for: iTime
1628           } // for: iCol
1629         } // for: iRow
1630
1631       } // if:  doMerge
1632
1633       mergeSDigitsManager->RemoveDigits(iDet);
1634       mergeSDigitsManager->RemoveDictionaries(iDet);
1635   
1636       if (fCompress) {
1637         digitsA->Compress(0); 
1638         for (iDict = 0; iDict < kNDict; iDict++) {                                     
1639           dictionaryA[iDict]->Compress();
1640         }
1641       }
1642       
1643     } // for: detectors    
1644       
1645     // The next set of s-digits
1646     mergeSDigitsManager = (AliTRDdigitsManager *) fSDigitsManagerList->After(mergeSDigitsManager);
1647     
1648   } // while: mergeDigitsManagers
1649   
1650   return kTRUE;
1651
1652 }
1653
1654 //_____________________________________________________________________________
1655 Bool_t AliTRDdigitizer::ConvertSDigits()
1656 {
1657   //
1658   // Converts s-digits to normal digits
1659   //
1660
1661   AliTRDarraySignal *digitsIn = 0x0;
1662
1663   if (!fSDigitsManager->HasSDigits()) {
1664     AliError("No s-digits in digits manager");
1665     return kFALSE;
1666   }
1667
1668   // Loop through the detectors
1669   for (Int_t det = 0; det < AliTRDgeometry::Ndet(); det++) {
1670
1671     // Get the merged s-digits (signals)
1672     digitsIn = (AliTRDarraySignal *) fSDigitsManager->GetSDigits(det);
1673     if (!digitsIn->HasData()) {
1674       AliDebug(2,Form("No digits for det=%d",det));
1675       continue;
1676     }
1677     
1678     // Convert the merged sdigits to digits
1679     if (!Signal2ADC(det,digitsIn)) {
1680       continue;
1681     }
1682
1683     // Copy the dictionary information to the output array
1684     if (!CopyDictionary(det)) {
1685       continue;
1686     }
1687
1688     // Delete 
1689     fSDigitsManager->RemoveDigits(det);
1690     fSDigitsManager->RemoveDictionaries(det);
1691
1692     // Run digital processing
1693     RunDigitalProcessing(det);
1694
1695     // Compress the arrays
1696     CompressOutputArrays(det);
1697
1698   } // for: detector numbers
1699
1700
1701   return kTRUE;
1702
1703 }
1704
1705 //_____________________________________________________________________________
1706 Bool_t AliTRDdigitizer::CopyDictionary(Int_t det)
1707 {
1708   //
1709   // Copies the dictionary information from the s-digits arrays
1710   // to the output arrays
1711   //
1712
1713   AliTRDcalibDB     *calibration = AliTRDcalibDB::Instance();
1714   if (!calibration) {
1715     AliFatal("Could not get calibration object");
1716     return kFALSE;
1717   }
1718
1719   AliDebug(1,Form("Start copying dictionaries for detector=%d",det));
1720
1721   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1722   AliTRDarrayDictionary *dictionaryIn[kNDict];
1723   AliTRDarrayDictionary *dictionaryOut[kNDict];
1724
1725   Int_t nRowMax    = fGeo->GetPadPlane(det)->GetNrows();
1726   Int_t nColMax    = fGeo->GetPadPlane(det)->GetNcols();
1727   Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
1728
1729   Int_t row  = 0;
1730   Int_t col  = 0;
1731   Int_t time = 0;
1732   Int_t dict = 0;
1733
1734   for (dict = 0; dict < kNDict; dict++) {
1735
1736     dictionaryIn[dict]  = (AliTRDarrayDictionary *) fSDigitsManager->GetDictionary(det,dict);
1737     dictionaryIn[dict]->Expand();
1738     dictionaryOut[dict] = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
1739     dictionaryOut[dict]->Allocate(nRowMax,nColMax,nTimeTotal);
1740
1741     for (row = 0; row < nRowMax; row++) {
1742       for (col = 0; col < nColMax; col++) {
1743         for (time = 0; time < nTimeTotal; time++) {
1744           Int_t track = dictionaryIn[dict]->GetData(row,col,time);
1745           dictionaryOut[dict]->SetData(row,col,time,track);
1746         } // for: time
1747       } // for: col
1748     } // for: row
1749     
1750   } // for: dictionaries
1751   
1752   return kTRUE;
1753
1754 }
1755
1756 //_____________________________________________________________________________
1757 void AliTRDdigitizer::CompressOutputArrays(Int_t det)
1758 {
1759   //
1760   // Compress the output arrays
1761   //
1762
1763   const Int_t kNDict = AliTRDdigitsManager::kNDict;
1764   AliTRDarrayDictionary *dictionary = 0x0;
1765
1766   if (fCompress) {
1767
1768     if (!fSDigits) {
1769       AliTRDarrayADC *digits = 0x0;  
1770       digits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det);
1771       digits->Compress();
1772     }
1773
1774     if (fSDigits) {
1775       AliTRDarraySignal *digits = 0x0; 
1776       digits = (AliTRDarraySignal *) fDigitsManager->GetSDigits(det);
1777       digits->Compress(0);
1778     }
1779
1780     for (Int_t dict = 0; dict < kNDict; dict++) {
1781       dictionary = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(det,dict);
1782       dictionary->Compress();
1783     }
1784
1785   }
1786
1787 }
1788
1789 //_____________________________________________________________________________
1790 Bool_t AliTRDdigitizer::WriteDigits() const
1791 {
1792   //
1793   // Writes out the TRD-digits and the dictionaries
1794   //
1795
1796   // Write parameters
1797   fRunLoader->CdGAFile();
1798
1799   // Store the digits and the dictionary in the tree
1800   return fDigitsManager->WriteDigits();
1801
1802 }
1803
1804 //_____________________________________________________________________________
1805 void AliTRDdigitizer::InitOutput(Int_t iEvent)
1806 {
1807   //
1808   // Initializes the output branches
1809   //
1810
1811   fEvent = iEvent;
1812    
1813   if (!fRunLoader) {
1814     AliError("Run Loader is NULL");
1815     return;  
1816   }
1817
1818   AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
1819   if (!loader) {
1820     AliError("Can not get TRD loader from Run Loader");
1821     return;
1822   }
1823
1824   TTree *tree = 0;
1825   
1826   if (fSDigits) { 
1827     // If we produce SDigits
1828     tree = loader->TreeS();
1829     if (!tree) {
1830       loader->MakeTree("S");
1831       tree = loader->TreeS();
1832     }
1833   }
1834   else {
1835     // If we produce Digits
1836     tree = loader->TreeD();
1837     if (!tree) {
1838       loader->MakeTree("D");
1839       tree = loader->TreeD();
1840     }
1841   }
1842   fDigitsManager->SetEvent(iEvent);
1843   fDigitsManager->MakeBranch(tree);
1844
1845 }
1846   
1847 //_____________________________________________________________________________
1848 Int_t AliTRDdigitizer::Diffusion(Float_t vdrift, Double_t absdriftlength
1849                                , Double_t &lRow, Double_t &lCol, Double_t &lTime)
1850 {
1851   //
1852   // Applies the diffusion smearing to the position of a single electron.
1853   // Depends on absolute drift length.
1854   //
1855   
1856   Float_t diffL = 0.0;
1857   Float_t diffT = 0.0;
1858
1859   if (AliTRDCommonParam::Instance()->GetDiffCoeff(diffL,diffT,vdrift)) {
1860
1861     Float_t driftSqrt = TMath::Sqrt(absdriftlength);
1862     Float_t sigmaT    = driftSqrt * diffT;
1863     Float_t sigmaL    = driftSqrt * diffL;
1864     lRow  = gRandom->Gaus(lRow ,sigmaT);
1865     lCol  = gRandom->Gaus(lCol ,sigmaT * GetLorentzFactor(vdrift));
1866     lTime = gRandom->Gaus(lTime,sigmaL * GetLorentzFactor(vdrift));
1867
1868     return 1;
1869
1870   }
1871   else {
1872
1873     return 0;
1874
1875   }
1876
1877 }
1878
1879 //_____________________________________________________________________________
1880 Float_t AliTRDdigitizer::GetLorentzFactor(Float_t vd)
1881 {
1882   //
1883   // Returns the Lorentz factor
1884   //
1885
1886   Double_t omegaTau      = AliTRDCommonParam::Instance()->GetOmegaTau(vd);
1887   Double_t lorentzFactor = 1.0;
1888   if (AliTRDCommonParam::Instance()->ExBOn()) {
1889     lorentzFactor = 1.0 / (1.0 + omegaTau*omegaTau);
1890   }
1891
1892   return lorentzFactor;
1893
1894 }
1895   
1896 //_____________________________________________________________________________
1897 Int_t AliTRDdigitizer::ExB(Float_t vdrift, Double_t driftlength, Double_t &lCol)
1898 {
1899   //
1900   // Applies E x B effects to the position of a single electron.
1901   // Depends on signed drift length.
1902   //
1903
1904   lCol = lCol 
1905        + AliTRDCommonParam::Instance()->GetOmegaTau(vdrift) 
1906        * driftlength;
1907
1908   return 1;
1909
1910 }
1911
1912 //_____________________________________________________________________________
1913 void AliTRDdigitizer::RunDigitalProcessing(Int_t det)
1914 {
1915   //
1916   // Run the digital processing in the TRAP
1917   //
1918
1919   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
1920
1921   //Create and initialize the mcm object 
1922   AliTRDmcmSim* mcmfast = new AliTRDmcmSim(); 
1923
1924   AliTRDarrayADC *digits = fDigitsManager->GetDigits(det);
1925   if (!digits)
1926     return;
1927
1928   //Call the methods in the mcm class using the temporary array as input  
1929   for(Int_t rob = 0; rob < digits->GetNrow() / 2; rob++)
1930   {
1931     for(Int_t mcm = 0; mcm < 16; mcm++)
1932     {
1933       mcmfast->Init(det, rob, mcm); 
1934       mcmfast->SetData(digits, fDigitsManager);
1935       mcmfast->Filter();
1936       if (feeParam->GetTracklet()) {
1937         mcmfast->Tracklet();
1938         mcmfast->StoreTracklets();
1939       }
1940       mcmfast->ZSMapping();
1941       mcmfast->WriteData(digits);
1942     }
1943   }
1944
1945   delete mcmfast;
1946
1947 }