Coverity fix.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSDigitizer.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 /* History of cvs commits:
19  *
20  * $Log: AliPHOSDigitizer.cxx,v $
21  * Revision 1.104  2007/12/18 09:08:18  hristov
22  * Splitting of the QA maker into simulation and reconstruction dependent parts (Yves)
23  *
24  * Revision 1.103  2007/11/07 11:25:06  schutz
25  * Comment out the QA checking before starting digitization
26  *
27  * Revision 1.102  2007/10/19 18:04:29  schutz
28  * The standalone QA data maker is called from AliSimulation and AliReconstruction outside the event loop; i.e. re-reading the data. The QA data making in the event loop has been commented out.
29  *
30  * Revision 1.101  2007/10/14 21:08:10  schutz
31  * Introduced the checking of QA results from previous step before entering the event loop
32  *
33  * Revision 1.100  2007/10/10 09:05:10  schutz
34  * Changing name QualAss to QA
35  *
36  * Revision 1.99  2007/09/30 17:08:20  schutz
37  * Introducing the notion of QA data acquisition cycle (needed by online)
38  *
39  * Revision 1.98  2007/09/26 14:22:17  cvetan
40  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
41  *
42  * Revision 1.97  2007/08/07 14:12:03  kharlov
43  * Quality assurance added (Yves Schutz)
44  *
45  * Revision 1.96  2007/04/28 10:43:36  policheh
46  * Dead channels simulation: digit energy sets to 0.
47  *
48  * Revision 1.95  2007/04/10 07:20:52  kharlov
49  * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
50  *
51  * Revision 1.94  2007/02/01 10:34:47  hristov
52  * Removing warnings on Solaris x86
53  *
54  * Revision 1.93  2006/10/17 13:17:01  kharlov
55  * Replace AliInfo by AliDebug
56  *
57  * Revision 1.92  2006/08/28 10:01:56  kharlov
58  * Effective C++ warnings fixed (Timur Pocheptsov)
59  *
60  * Revision 1.91  2006/04/29 20:25:30  hristov
61  * Decalibration is implemented (Yu.Kharlov)
62  *
63  * Revision 1.90  2006/04/22 10:30:17  hristov
64  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
65  *
66  * Revision 1.89  2006/04/11 15:22:59  hristov
67  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
68  *
69  * Revision 1.88  2006/03/13 14:05:43  kharlov
70  * Calibration objects for EMC and CPV
71  *
72  * Revision 1.87  2005/08/24 15:33:49  kharlov
73  * Calibration data for raw digits
74  *
75  * Revision 1.86  2005/07/12 20:07:35  hristov
76  * Changes needed to run simulation and reconstrruction in the same AliRoot session
77  *
78  * Revision 1.85  2005/05/28 14:19:04  schutz
79  * Compilation warnings fixed by T.P.
80  *
81  */
82
83 //_________________________________________________________________________
84 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute) 
85 //////////////////////////////////////////////////////////////////////////////
86 // This class performs digitization of Summable digits (in the PHOS case it is just
87 // the sum of contributions from all primary particles into a given cell). 
88 // In addition it performs mixing of summable digits from different events.
89 // The name of the class is also the title of the branch that will contain 
90 // the created SDigits
91 // The title of the class is the name of the file that contains the hits from
92 // which the SDigits are created
93 //
94 // For each event two branches are created in TreeD:
95 //   "PHOS" - list of digits
96 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
97 //
98 // Note, that one can set a title for new digits branch, and repeat digitization with
99 // another set of parameters.
100 //
101 // Use case:
102 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
103 // root[1] d->Digitize()             
104 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
105 //                       //Digitizes SDigitis in all events found in file galice.root 
106 //
107 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;  
108 //                       // Will read sdigits from galice1.root
109 // root[3] d1->MixWith("galice2.root")       
110 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
111 //                       // Reads another set of sdigits from galice2.root
112 // root[3] d1->MixWith("galice3.root")       
113 //                       // Reads another set of sdigits from galice3.root
114 // root[4] d->Digitize("deb timing")    
115 //                       // Reads SDigits from files galice1.root, galice2.root ....
116 //                       // mixes them and stores produced Digits in file galice1.root          
117 //                       // deb - prints number of produced digits
118 //                       // deb all - prints list of produced digits
119 //                       // timing  - prints time used for digitization
120 //
121
122 // --- ROOT system ---
123 #include "TTree.h"
124 #include "TSystem.h"
125 #include "TBenchmark.h"
126 #include "TRandom.h"
127 #include "TMath.h"
128
129 // --- Standard library ---
130
131 // --- AliRoot header files ---
132 #include <TGeoManager.h>                                                                                                                   
133 #include "AliLog.h"
134 #include "AliDigitizationInput.h"
135 #include "AliPHOSDigit.h"
136 #include "AliPHOSDigitizer.h"
137 #include "AliPHOSGeometry.h"
138 #include "AliPHOSTick.h"
139 #include "AliPHOSSimParam.h"
140 #include "AliPHOSCalibData.h"
141 #include "AliRunLoader.h"
142 #include "AliPHOSLoader.h"
143 #include "AliPHOSPulseGenerator.h"
144
145 ClassImp(AliPHOSDigitizer)
146
147
148 //____________________________________________________________________________ 
149 AliPHOSDigitizer::AliPHOSDigitizer() :
150   AliDigitizer("",""),
151   fDefaultInit(kTRUE),
152   fDigitsInRun(0),
153   fInit(kFALSE),
154   fInput(0),
155   fInputFileNames(0x0),
156   fEventNames(0x0),
157   fEmcCrystals(0),
158   fEventFolderName(""),
159   fFirstEvent(0),
160   fLastEvent(0), 
161   fcdb(0x0),
162   fEventCounter(0),
163   fPulse(0),
164   fADCValuesLG(0),
165   fADCValuesHG(0)
166 {
167   // ctor
168   InitParameters() ; 
169   fDigInput = 0 ;                     // We work in the standalong mode
170 }
171
172 //____________________________________________________________________________ 
173 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName, 
174                                    TString eventFolderName):
175   AliDigitizer("PHOSDigitizer", alirunFileName), 
176   fDefaultInit(kFALSE),
177   fDigitsInRun(0),
178   fInit(kFALSE),
179   fInput(0),
180   fInputFileNames(0x0),
181   fEventNames(0x0),
182   fEmcCrystals(0),
183   fEventFolderName(eventFolderName),
184   fFirstEvent(0),
185   fLastEvent(0), 
186   fcdb(0x0),
187   fEventCounter(0),
188   fPulse(0),
189   fADCValuesLG(0),
190   fADCValuesHG(0)
191 {
192   // ctor
193   InitParameters() ; 
194   Init() ;
195   fDefaultInit = kFALSE ; 
196   fDigInput = 0 ;                     // We work in the standalone mode
197   fcdb = new AliPHOSCalibData(-1);
198 }
199
200 //____________________________________________________________________________ 
201 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) : 
202   AliDigitizer(d),
203   fDefaultInit(d.fDefaultInit),
204   fDigitsInRun(d.fDigitsInRun),
205   fInit(d.fInit),
206   fInput(d.fInput),
207   fInputFileNames(0x0),
208   fEventNames(0x0),
209   fEmcCrystals(d.fEmcCrystals),
210   fEventFolderName(d.fEventFolderName),
211   fFirstEvent(d.fFirstEvent),
212   fLastEvent(d.fLastEvent), 
213   fcdb (new AliPHOSCalibData(-1)), 
214   fEventCounter(0),
215   fPulse(new AliPHOSPulseGenerator),
216   fADCValuesLG(0),
217   fADCValuesHG(0)
218 {
219   // copyy ctor 
220   SetName(d.GetName()) ; 
221   SetTitle(d.GetTitle()) ; 
222   for (Int_t iInput=0; iInput<fInput; iInput++) {
223     fInputFileNames[iInput] = d.fInputFileNames[iInput];
224     fEventNames[iInput]     = d.fEventNames[iInput];
225   }
226   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
227   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
228 }
229
230 //____________________________________________________________________________ 
231 AliPHOSDigitizer::AliPHOSDigitizer(AliDigitizationInput * rd) :
232   AliDigitizer(rd,"PHOSDigitizer"),
233   fDefaultInit(kFALSE),
234   fDigitsInRun(0),
235   fInit(kFALSE),
236   fInput(0),
237   fInputFileNames(0x0),
238   fEventNames(0x0),
239   fEmcCrystals(0),
240   fEventFolderName(fDigInput->GetInputFolderName(0)),
241   fFirstEvent(0),
242   fLastEvent(0), 
243   fcdb (0x0), 
244   fEventCounter(0),
245   fPulse(0),
246   fADCValuesLG(0),
247   fADCValuesHG(0)
248
249 {
250   // ctor Init() is called by RunDigitizer
251   fDigInput = rd ; 
252   SetTitle(static_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
253   InitParameters() ; 
254   fDefaultInit = kFALSE ; 
255   fcdb = new AliPHOSCalibData(-1);
256 }
257
258 //____________________________________________________________________________ 
259   AliPHOSDigitizer::~AliPHOSDigitizer()
260 {
261   // dtor  
262   delete [] fInputFileNames ; 
263   delete [] fEventNames ; 
264
265   delete fPulse;
266   delete [] fADCValuesLG;
267   delete [] fADCValuesHG;
268
269   if(fcdb){ delete fcdb ; fcdb=0;} 
270
271 }
272
273 //____________________________________________________________________________
274 void AliPHOSDigitizer::Digitize(Int_t event) 
275
276   
277   // Makes the digitization of the collected summable digits.
278   //  It first creates the array of all PHOS modules
279   //  filled with noise (different for EMC, and CPV) and
280   //  then adds contributions from SDigits. 
281   // This design avoids scanning over the list of digits to add 
282   // contribution to new SDigits only.
283
284   Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
285
286   //First stream 
287   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
288   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
289
290   Int_t readEvent = event ; 
291   if (fDigInput) 
292     readEvent = static_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ; 
293   AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
294                   readEvent, GetTitle(), fEventFolderName.Data())) ; 
295   rl->GetEvent(readEvent) ;
296   phosLoader->CleanSDigits() ; 
297   phosLoader->LoadSDigits("READ") ;
298
299   //Prepare Output
300   TClonesArray * digits = phosLoader->Digits() ;
301   if( !digits ) {
302     phosLoader->MakeDigitsArray() ; 
303     digits = phosLoader->Digits() ;
304   }
305  
306   digits->Clear() ;
307
308   //
309   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
310   //Making digits with noise, first EMC
311   //Check which PHOS modules are present
312   Bool_t isPresent[5] ;
313   TString volpath ;
314   Int_t nmod=0 ;
315   for(Int_t i=0; i<5; i++){
316     volpath = "/ALIC_1/PHOS_";
317     volpath += i+1;
318     if (gGeoManager->CheckPath(volpath.Data())) {
319       isPresent[i]=1 ;
320       nmod++ ;
321     }
322     else{
323       isPresent[i]=0 ;
324     }
325   }
326
327   Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
328   
329   Int_t nCPV ;
330   Int_t absID ;
331   
332   //check if CPV exists
333   Bool_t isCPVpresent=0 ;
334   for(Int_t i=1; i<=5 && !isCPVpresent; i++){
335     volpath = "/ALIC_1/PHOS_";
336     volpath += i;
337     volpath += "/PCPV_1";
338     if (gGeoManager->CheckPath(volpath.Data())) 
339       isCPVpresent=1 ;
340   } 
341   
342   if(isCPVpresent){
343     nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
344   }
345   else{
346      nCPV = nEMC ;
347   }  
348
349   digits->Expand(nCPV) ;
350
351   //take all the inputs to add together and load the SDigits
352   TObjArray * sdigArray = new TObjArray(fInput) ;
353   sdigArray->AddAt(phosLoader->SDigits(), 0) ;
354  
355   for(Int_t i = 1 ; i < fInput ; i++){
356     TString tempo(fEventNames[i]) ; 
357     tempo += i ;
358     AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
359     if(!rl2){ 
360       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
361       if (!rl2) {
362         Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
363         return ;
364       }
365       rl2->LoadHeader();
366     }
367     AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
368  
369     if(fDigInput){ 
370       readEvent = static_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ; 
371     }
372     TClonesArray * digs ;
373     if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
374       AliInfo(Form("Adding event %d from input stream %d %s %s", 
375                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ; 
376        rl2->GetEvent(readEvent) ;
377        phosLoader2->CleanDigits() ;
378        phosLoader2->LoadDigits("READ") ;
379        digs = phosLoader2->Digits() ;
380        toMakeNoise=0 ; //Do not add noise, it is already in stream
381     }
382     else{
383       AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
384                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
385        rl2->GetEvent(readEvent) ;
386        phosLoader2->CleanSDigits() ; 
387        phosLoader2->LoadSDigits("READ") ;
388        digs = phosLoader2->SDigits() ;
389     } 
390     sdigArray->AddAt(digs, i) ;
391   }
392
393   //Find the first crystal with signal
394   Int_t nextSig = 200000 ; 
395   TClonesArray * sdigits ;  
396   for(Int_t i = 0 ; i < fInput ; i++){
397     sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
398     if ( !sdigits->GetEntriesFast() )
399       continue ; 
400     Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
401     if(curNext < nextSig) 
402       nextSig = curNext ;
403   }
404   
405   TArrayI index(fInput) ;
406   index.Reset() ;  //Set all indexes to zero
407   
408   AliPHOSDigit * digit ;
409   AliPHOSDigit * curSDigit ;
410   
411 //  TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
412   
413   //Put Noise contribution
414   Double_t apdNoise = 0. ;
415   if(toMakeNoise)
416      apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ; 
417
418   Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
419   Int_t idigit= 0;
420   for(Int_t imod=0; imod<5; imod++){
421     if(!isPresent[imod])
422       continue ;
423     Int_t firstAbsId=imod*emcpermod+1 ;
424     Int_t lastAbsId =(imod+1)*emcpermod ; 
425     for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
426       Float_t noise = gRandom->Gaus(0.,apdNoise) ; 
427       new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
428       //look if we have to add signal?
429       digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
430       idigit++ ;
431     
432       if(absID==nextSig){
433       //Add SDigits from all inputs 
434 //      ticks->Clear() ;
435 //      Int_t contrib = 0 ;
436
437 //New Timing model is necessary
438 //      Float_t a = digit->GetEnergy() ;
439 //      Float_t b = TMath::Abs( a / fTimeSignalLength) ;
440 //      //Mark the beginning of the signal
441 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);  
442 //      //Mark the end of the signal     
443 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b); 
444
445 // Calculate time as time of the largest digit
446         Float_t time = digit->GetTime() ;
447         Float_t eTime= digit->GetEnergy() ;
448       
449         //loop over inputs
450         for(Int_t i = 0 ; i < fInput ; i++){
451           if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
452             curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;       
453             if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
454               curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
455               curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
456             }
457           }
458           else
459             curSDigit = 0 ;
460           //May be several digits will contribute from the same input
461           while(curSDigit && curSDigit->GetId() == absID){         
462             //Shift primary to separate primaries belonging different inputs
463             Int_t primaryoffset ;
464             if(fDigInput)
465               primaryoffset = fDigInput->GetMask(i) ; 
466             else
467               primaryoffset = 10000000*i ;
468             curSDigit->ShiftPrimary(primaryoffset) ;
469           
470 //New Timing model is necessary
471 //        a = curSDigit->GetEnergy() ;
472 //        b = a /fTimeSignalLength ;
473 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);  
474 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
475           if(curSDigit->GetEnergy()>eTime){
476             eTime=curSDigit->GetEnergy() ;
477             time=curSDigit->GetTime() ;
478           }
479             *digit += *curSDigit ;  //add energies
480
481             index[i]++ ;
482             if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
483               curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;     
484             else
485               curSDigit = 0 ;
486           }
487         }
488       
489         //calculate and set time
490 //New Timing model is necessary
491 //      Float_t time = FrontEdgeTime(ticks) ;
492         digit->SetTime(time) ;
493       
494         //Find next signal module
495         nextSig = 200000 ;
496         for(Int_t i = 0 ; i < fInput ; i++){
497           sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
498           Int_t curNext = nextSig ;
499           if(sdigits->GetEntriesFast() > index[i] ){
500             curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
501           }
502           if(curNext < nextSig) nextSig = curNext ;
503         }
504       }
505     }
506   }
507
508
509   //Apply non-linearity
510   if(AliPHOSSimParam::GetInstance()->IsCellNonlinearityOn()){ //Apply non-lineairyt on cell level
511     const Double_t aNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyA() ;
512     const Double_t bNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyB() ;
513     const Double_t cNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyC() ;
514     for(Int_t i = 0 ; i < nEMC ; i++){
515       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
516       Double_t e= digit->GetEnergy() ;
517       // version(1)      digit->SetEnergy(e*(1+a*TMath::Exp(-e/b))) ;
518       digit->SetEnergy(e*cNL*(1.+aNL*TMath::Exp(-e*e/2./bNL/bNL))) ; //Better agreement with data...
519     }  
520   }
521
522
523   //distretize energy if necessary
524   if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
525     Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
526     for(Int_t i = 0 ; i < nEMC ; i++){                                                                                                       
527       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
528       digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
529     } 
530   }
531  
532   //Apply decalibration if necessary
533   for(Int_t i = 0 ; i < nEMC ; i++){
534     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
535     Decalibrate(digit) ;
536   }
537   
538 //  ticks->Delete() ;
539 //  delete ticks ;
540   
541   //Now CPV digits (different noise and no timing)
542   Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
543   Int_t nEMCtotal=emcpermod*5 ;
544   Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
545   if(isCPVpresent){  //CPV is present in current geometry
546     for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
547       if(!isPresent[imod])
548         continue ;
549       Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
550       Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
551       for(absID = firstAbsId; absID <= lastAbsId; absID++){
552         Float_t noise = gRandom->Gaus(0., cpvNoise) ; 
553         new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
554         idigit++ ;
555         //look if we have to add signal?
556         if(absID==nextSig){
557           digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
558           //Add SDigits from all inputs
559           for(Int_t i = 0 ; i < fInput ; i++){
560              if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
561                curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;   
562              else
563                curSDigit = 0 ;
564
565              //May be several digits will contribute from the same input
566              while(curSDigit && curSDigit->GetId() == absID){      
567                //Shift primary to separate primaries belonging different inputs
568                Int_t primaryoffset ;
569                if(fDigInput)
570                  primaryoffset = fDigInput->GetMask(i) ; 
571                else
572                  primaryoffset = 10000000*i ;
573                curSDigit->ShiftPrimary(primaryoffset) ;
574
575                //add energies
576                *digit += *curSDigit ;  
577                index[i]++ ;
578                if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
579                  curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;        
580                else
581                  curSDigit = 0 ;
582              }
583           }
584
585           //Find next signal module
586           nextSig = 200000 ;
587           for(Int_t i = 0 ; i < fInput ; i++){
588             sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
589             Int_t curNext = nextSig ;
590             if(sdigits->GetEntriesFast() > index[i] )
591               curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
592             if(curNext < nextSig) nextSig = curNext ;
593           }
594       
595         }
596       }
597     }
598   }
599
600   delete sdigArray ; //We should not delete its contents 
601   
602   Int_t relId[4];
603
604   //set amplitudes in bad channels to zero
605
606   for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
607     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
608     geom->AbsToRelNumbering(digit->GetId(),relId);
609     if(relId[1] == 0) // Emc
610       if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.); 
611   }
612
613   //remove digits below thresholds
614   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
615   for(Int_t i = 0 ; i < nEMC ; i++){
616     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
617
618     if(digit->GetEnergy() < emcThreshold){
619       digits->RemoveAt(i) ;
620       continue ;
621     }
622
623     geom->AbsToRelNumbering(digit->GetId(),relId);
624
625 //    digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
626
627     Float_t tres = TimeResolution(digit->GetEnergy()) ; 
628     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
629
630     fPulse->Reset();
631     fPulse->SetAmplitude(digit->GetEnergy()/
632                          fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
633     fPulse->SetTZero(digit->GetTimeR());
634     fPulse->MakeSamples();
635     fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ; 
636     Int_t nSamples = fPulse->GetRawFormatTimeBins();
637     digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
638     digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
639   }
640
641   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
642   for(Int_t i = nEMC; i < nCPV ; i++){
643     if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
644       digits->RemoveAt(i) ;
645   } 
646     
647   digits->Compress() ;  
648   Int_t ndigits = digits->GetEntriesFast() ;
649   
650   //Set indexes in list of digits and make true digitization of the energy
651   for (Int_t i = 0 ; i < ndigits ; i++) { 
652     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; 
653     digit->SetIndexInList(i) ;     
654     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
655       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
656     }
657   }
658
659 }
660 //____________________________________________________________________________
661 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
662   //Apply calibration
663   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
664
665   //Determine rel.position of the cell absolute ID
666   Int_t relId[4];
667   geom->AbsToRelNumbering(absId,relId);
668   Int_t module=relId[0];
669   Int_t row   =relId[2];
670   Int_t column=relId[3];
671   if(relId[1]==0){ //This Is EMC
672     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
673     return amp*calibration ;
674   }
675   return 0 ;
676 }
677 //____________________________________________________________________________
678 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
679 {
680   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
681
682   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
683
684   //Determine rel.position of the cell absolute ID
685   Int_t relId[4];
686   geom->AbsToRelNumbering(digit->GetId(),relId);
687   Int_t module=relId[0];
688   Int_t row   =relId[2];
689   Int_t column=relId[3];
690   if(relId[1]==0){ //This Is EMC
691     Float_t decalib     = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
692     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
693     Float_t energy = digit->GetEnergy()/calibration;
694     digit->SetEnergy(energy); //Now digit measures E in ADC counts
695     Float_t time = digit->GetTime() ;
696     time-=fcdb->GetTimeShiftEmc(module,column,row);
697     digit->SetTime(time) ;
698   }
699 }
700 //____________________________________________________________________________
701 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
702   //Apply time calibration
703   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
704
705   //Determine rel.position of the cell absolute ID
706   Int_t relId[4];
707   geom->AbsToRelNumbering(absId,relId);
708   Int_t module=relId[0];
709   Int_t row   =relId[2];
710   Int_t column=relId[3];
711   time += fcdb->GetTimeShiftEmc(module,column,row);
712   return time ;
713 }
714 //____________________________________________________________________________
715 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
716 {
717   // Returns digitized value of the CPV charge in a pad absId
718
719   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
720
721   //Determine rel.position of the cell absId
722   Int_t relId[4];
723   geom->AbsToRelNumbering(absId,relId);
724   Int_t module=relId[0];
725   Int_t row   =relId[2];
726   Int_t column=relId[3];
727   
728   Int_t channel = 0;
729   
730   if(absId > fEmcCrystals){ //digitize CPV only
731
732     //reading calibration data for cell absId.
733     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
734     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
735
736     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;       
737     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
738     if(channel > nMax ) channel = nMax ;
739   }
740   return channel ;
741 }
742
743 //____________________________________________________________________________
744 void AliPHOSDigitizer::Digitize(Option_t *option) 
745
746   // Steering method to process digitization for events
747   // in the range from fFirstEvent to fLastEvent.
748   // This range is optionally set by SetEventRange().
749   // if fLastEvent=-1, then process events until the end.
750   // by default fLastEvent = fFirstEvent (process only one event)
751
752   if (!fInit) { // to prevent overwrite existing file
753     AliError(Form("Give a version name different from %s", 
754                   fEventFolderName.Data() )) ;
755     return ;
756   }   
757
758   if (strstr(option,"print")) {
759     Print();
760     return ; 
761   }
762   
763  if(strstr(option,"tim"))
764      gBenchmark->Start("PHOSDigitizer");
765   
766   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
767   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
768   
769   if (fLastEvent == -1) 
770     fLastEvent = rl->GetNumberOfEvents() - 1 ;
771   else if (fDigInput) 
772     fLastEvent = fFirstEvent ; 
773  
774   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
775   
776   Int_t ievent ;
777
778   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
779     fEventCounter++ ; 
780
781     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
782
783     WriteDigits() ;
784
785     if(strstr(option,"deb"))
786       PrintDigits(option);
787     
788     //increment the total number of Digits per run 
789     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;  
790  }
791  
792   if(strstr(option,"tim")){
793     gBenchmark->Stop("PHOSDigitizer");
794     TString message ; 
795     message = "  took %f seconds for Digitizing %f seconds per event\n" ; 
796     AliInfo(Form( message.Data(), 
797          gBenchmark->GetCpuTime("PHOSDigitizer"), 
798          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents )); 
799   } 
800 }
801 //____________________________________________________________________________ 
802 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
803   //calculate TOF resolution using beam-test resutls
804   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
805   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
806   return TMath::Sqrt(a*a+b*b/e) ;
807 }
808
809 ////____________________________________________________________________________ 
810 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
811 //{
812 //  // Returns the shortest time among all time ticks
813 //
814 //  ticks->Sort() ; //Sort in accordance with times of ticks
815 //  TIter it(ticks) ;
816 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
817 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
818 //
819 //  AliPHOSTick * t ;  
820 //  while((t=(AliPHOSTick*) it.Next())){
821 //    if(t->GetTime() < time)  //This tick starts before crossing
822 //      *ctick+=*t ;
823 //    else
824 //      return time ;
825 //
826 //    time = ctick->CrossingTime(fTimeThreshold) ;    
827 //  }
828 //  return time ;
829 //}
830
831 //____________________________________________________________________________ 
832 Bool_t AliPHOSDigitizer::Init()
833 {
834   // Makes all memory allocations
835   fInit = kTRUE ; 
836
837   AliPHOSGeometry *geom;
838   if (!(geom = AliPHOSGeometry::GetInstance())) 
839         geom = AliPHOSGeometry::GetInstance("IHEP","");
840 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
841
842   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
843   
844   fFirstEvent = 0 ; 
845   fLastEvent = fFirstEvent ; 
846   if (fDigInput) 
847     fInput = fDigInput->GetNinputs() ; 
848   else 
849     fInput           = 1 ;
850   
851   fInputFileNames  = new TString[fInput] ; 
852   fEventNames      = new TString[fInput] ; 
853   fInputFileNames[0] = GetTitle() ; 
854   fEventNames[0]     = fEventFolderName.Data() ; 
855   Int_t index ; 
856   for (index = 1 ; index < fInput ; index++) {
857     fInputFileNames[index] = static_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0); 
858     TString tempo = fDigInput->GetInputFolderName(index) ;
859     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fDigInput
860   }
861
862   //to prevent cleaning of this object while GetEvent is called
863   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
864   if(!rl){
865     AliRunLoader::Open(GetTitle(), fEventFolderName) ; 
866   }
867   return fInit ; 
868 }
869
870 //____________________________________________________________________________ 
871 void AliPHOSDigitizer::InitParameters()
872 {
873   // Set initial parameters Digitizer
874
875   fDigitsInRun  = 0 ; 
876   SetEventRange(0,-1) ;
877   fPulse = new AliPHOSPulseGenerator();
878   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
879   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
880     
881 }
882
883 //__________________________________________________________________
884 void AliPHOSDigitizer::Print(const Option_t *)const 
885 {
886   // Print Digitizer's parameters
887   AliInfo(Form("\n------------------- %s -------------", GetName() )) ; 
888   if( strcmp(fEventFolderName.Data(), "") != 0 ){
889     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
890     
891     Int_t nStreams ; 
892     if (fDigInput) 
893       nStreams =  GetNInputStreams() ;
894     else 
895       nStreams = fInput ; 
896     
897     Int_t index = 0 ;  
898     for (index = 0 ; index < nStreams ; index++) {  
899       TString tempo(fEventNames[index]) ; 
900       tempo += index ;
901       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ; 
902     }
903     
904  //   printf("\nWith following parameters:\n") ;
905  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ; 
906  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;  
907  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ; 
908  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;  
909     printf(" ---------------------------------------------------\n") ;   
910   }
911   else
912     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
913   
914 }
915
916 //__________________________________________________________________
917  void AliPHOSDigitizer::PrintDigits(Option_t * option)
918 {
919   // Print a table of digits
920   
921
922   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
923   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
924   TClonesArray * digits = phosLoader->Digits() ; 
925   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
926   
927   AliInfo(Form("%d", digits->GetEntriesFast())) ; 
928   printf("\nevent %d", gAlice->GetEvNumber()) ;
929   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;  
930
931  
932   if(strstr(option,"all")||strstr(option,"EMC")){  
933     //loop over digits
934     AliPHOSDigit * digit;
935     printf("\nEMC digits (with primaries):\n")  ;
936     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
937     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
938     Int_t index ;
939     for (index = 0 ; (index < digits->GetEntriesFast()) && 
940            (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
941       digit = (AliPHOSDigit * )  digits->At(index) ;
942       if(digit->GetNprimary() == 0) 
943         continue;
944 //       printf("%6d  %8d    %6.5e %4d      %2d :",
945 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
946       printf("%6d  %.4f    %6.5e %4d      %2d :",
947               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
948       Int_t iprimary;
949       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
950         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
951       }    
952       printf("\n") ;
953     }
954   }
955   
956   if(strstr(option,"all")||strstr(option,"CPV")){
957     
958     //loop over CPV digits
959     AliPHOSDigit * digit;
960     printf("\nCPV digits (with primaries):\n")  ;
961     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
962     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
963     Int_t index ;
964     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
965       digit = (AliPHOSDigit * )  digits->At(index) ;
966       if(digit->GetId() > maxEmc){
967         printf("%6d  %8d    %4d      %2d :",
968                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;  
969         Int_t iprimary;
970         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
971           printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
972         }    
973         printf("\n") ;
974       }
975     }
976   }
977
978 }
979
980 //__________________________________________________________________
981 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
982 {  // Calculates the time signal generated by noise
983   //PH  Info("TimeOfNoise", "Change me") ; 
984   return gRandom->Rndm() * 1.28E-5;
985 }
986
987 //__________________________________________________________________
988 void AliPHOSDigitizer::Unload() 
989 {  
990   
991   Int_t i ; 
992   for(i = 1 ; i < fInput ; i++){
993     TString tempo(fEventNames[i]) ; 
994     tempo += i ;
995     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ; 
996     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
997     phosLoader->UnloadSDigits() ; 
998   }
999   
1000   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1001   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1002   phosLoader->UnloadDigits() ; 
1003 }
1004
1005 //____________________________________________________________________________
1006 void AliPHOSDigitizer::WriteDigits()
1007 {
1008
1009   // Makes TreeD in the output file. 
1010   // Check if branch already exists: 
1011   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
1012   //      already existing branches. 
1013   //   else creates branch with Digits, named "PHOS", title "...",
1014   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
1015   //      and names of files, from which digits are made.
1016
1017   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1018   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1019  
1020   const TClonesArray * digits = phosLoader->Digits() ; 
1021   TTree * treeD = phosLoader->TreeD();
1022   if(!treeD){
1023     phosLoader->MakeTree("D");
1024     treeD = phosLoader->TreeD();
1025   }
1026   
1027   // -- create Digits branch
1028   Int_t bufferSize = 32000 ;    
1029   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1030   digitsBranch->SetTitle(fEventFolderName);
1031   digitsBranch->Fill() ;
1032   
1033   phosLoader->WriteDigits("OVERWRITE");
1034
1035   Unload() ; 
1036
1037 }