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