88dabf0aa5decac00834aaa8ff74126bb9ee60a2
[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 TTask 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 TTask is also the title of the branch that will contain 
90 // the created SDigits
91 // The title of the TTAsk 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->ExecuteTask()             
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->ExecuteTask("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 "AliRunDigitizer.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
144 ClassImp(AliPHOSDigitizer)
145
146
147 //____________________________________________________________________________ 
148 AliPHOSDigitizer::AliPHOSDigitizer() :
149   AliDigitizer("",""),
150   fDefaultInit(kTRUE),
151   fDigitsInRun(0),
152   fInit(kFALSE),
153   fInput(0),
154   fInputFileNames(0x0),
155   fEventNames(0x0),
156   fEmcCrystals(0),
157   fEventFolderName(""),
158   fFirstEvent(0),
159   fLastEvent(0), 
160   fcdb(0x0),
161   fEventCounter(0)
162 {
163   // ctor
164   InitParameters() ; 
165   fManager = 0 ;                     // We work in the standalong mode
166 }
167
168 //____________________________________________________________________________ 
169 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName, 
170                                    TString eventFolderName):
171   AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName), 
172   fDefaultInit(kFALSE),
173   fDigitsInRun(0),
174   fInit(kFALSE),
175   fInput(0),
176   fInputFileNames(0x0),
177   fEventNames(0x0),
178   fEmcCrystals(0),
179   fEventFolderName(eventFolderName),
180   fFirstEvent(0),
181   fLastEvent(0), 
182   fcdb(0x0),
183   fEventCounter(0)
184 {
185   // ctor
186   InitParameters() ; 
187   Init() ;
188   fDefaultInit = kFALSE ; 
189   fManager = 0 ;                     // We work in the standalong mode
190   fcdb = new AliPHOSCalibData(-1);
191 }
192
193 //____________________________________________________________________________ 
194 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) : 
195   AliDigitizer(d),
196   fDefaultInit(d.fDefaultInit),
197   fDigitsInRun(d.fDigitsInRun),
198   fInit(d.fInit),
199   fInput(d.fInput),
200   fInputFileNames(0x0),//?
201   fEventNames(0x0),//?
202   fEmcCrystals(d.fEmcCrystals),
203   fEventFolderName(d.fEventFolderName),
204   fFirstEvent(d.fFirstEvent),
205   fLastEvent(d.fLastEvent), 
206   fcdb (0x0), 
207   fEventCounter(0)
208 {
209   // copyy ctor 
210   SetName(d.GetName()) ; 
211   SetTitle(d.GetTitle()) ; 
212   fcdb = new AliPHOSCalibData(-1);
213 }
214
215 //____________________________________________________________________________ 
216 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
217   AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
218   fDefaultInit(kFALSE),
219   fDigitsInRun(0),
220   fInit(kFALSE),
221   fInput(0),
222   fInputFileNames(0x0),
223   fEventNames(0x0),
224   fEmcCrystals(0),
225   fEventFolderName(fManager->GetInputFolderName(0)),
226   fFirstEvent(0),
227   fLastEvent(0), 
228   fcdb (0x0), 
229   fEventCounter(0)
230
231 {
232   // ctor Init() is called by RunDigitizer
233   fManager = rd ; 
234   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
235   InitParameters() ; 
236   fDefaultInit = kFALSE ; 
237   fcdb = new AliPHOSCalibData(-1);
238 }
239
240 //____________________________________________________________________________ 
241   AliPHOSDigitizer::~AliPHOSDigitizer()
242 {
243   // dtor
244   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
245   if(rl){                                                                                                                               
246     AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));                                                                 
247     phosLoader->CleanDigitizer() ;
248   }
249
250   delete [] fInputFileNames ; 
251   delete [] fEventNames ; 
252
253   if(fcdb){ delete fcdb ; fcdb=0;} 
254
255 }
256
257 //____________________________________________________________________________
258 void AliPHOSDigitizer::Digitize(Int_t event) 
259
260   
261   // Makes the digitization of the collected summable digits.
262   //  It first creates the array of all PHOS modules
263   //  filled with noise (different for EMC, and CPV) and
264   //  then adds contributions from SDigits. 
265   // This design avoids scanning over the list of digits to add 
266   // contribution to new SDigits only.
267
268   Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
269
270   //First stream 
271   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
272   AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));                                                                
273
274   Int_t readEvent = event ; 
275   if (fManager) 
276     readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ; 
277   AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
278                   readEvent, GetTitle(), fEventFolderName.Data())) ; 
279   rl->GetEvent(readEvent) ;
280   phosLoader->CleanSDigits() ; 
281   phosLoader->LoadSDigits("READ") ;
282
283   //Prepare Output
284   TClonesArray * digits = phosLoader->Digits() ;
285   if( !digits ) {
286     phosLoader->MakeDigitsArray() ; 
287     digits = phosLoader->Digits() ;
288   }
289  
290   digits->Clear() ;
291
292   //
293   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
294   //Making digits with noise, first EMC
295   //Check which PHOS modules are present
296   Bool_t isPresent[5] ;
297   TString volpath ;
298   Int_t nmod=0 ;
299   for(Int_t i=0; i<5; i++){
300     volpath = "/ALIC_1/PHOS_";
301     volpath += i+1;
302     if (gGeoManager->CheckPath(volpath.Data())) {
303       isPresent[i]=1 ;
304       nmod++ ;
305     }
306     else{
307       isPresent[i]=0 ;
308     }
309   }
310
311   Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
312   
313   Int_t nCPV ;
314   Int_t absID ;
315   
316   //check if CPV exists
317   Bool_t isCPVpresent=0 ;
318   for(Int_t i=1; i<=5 && !isCPVpresent; i++){
319     volpath = "/ALIC_1/PHOS_";
320     volpath += i;
321     volpath += "/PCPV_1";
322     if (gGeoManager->CheckPath(volpath.Data())) 
323       isCPVpresent=1 ;
324   } 
325   
326   if(isCPVpresent){
327     nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
328   }
329   else{
330      nCPV = nEMC ;
331   }  
332
333   digits->Expand(nCPV) ;
334
335   //take all the inputs to add together and load the SDigits
336   TObjArray * sdigArray = new TObjArray(fInput) ;
337   sdigArray->AddAt(phosLoader->SDigits(), 0) ;
338  
339   for(Int_t i = 1 ; i < fInput ; i++){
340     TString tempo(fEventNames[i]) ; 
341     tempo += i ;
342     AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
343     if(!rl2){ 
344       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
345       if (!rl2) {
346         Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
347         return ;
348       }
349       rl2->LoadHeader();
350     }
351     AliPHOSLoader * phosLoader2 = dynamic_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
352  
353     if(fManager){ 
354       readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ; 
355     }
356     TClonesArray * digs ;
357     if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
358       AliInfo(Form("Adding event %d from input stream %d %s %s", 
359                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ; 
360        rl2->GetEvent(readEvent) ;
361        phosLoader2->CleanDigits() ;
362        phosLoader2->LoadDigits("READ") ;
363        digs = phosLoader2->Digits() ;
364        toMakeNoise=0 ; //Do not add noise, it is already in stream
365     }
366     else{
367       AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
368                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
369        rl2->GetEvent(readEvent) ;
370        phosLoader2->CleanSDigits() ; 
371        phosLoader2->LoadSDigits("READ") ;
372        digs = phosLoader2->SDigits() ;
373     } 
374     sdigArray->AddAt(digs, i) ;
375   }
376
377   //Find the first crystal with signal
378   Int_t nextSig = 200000 ; 
379   TClonesArray * sdigits ;  
380   for(Int_t i = 0 ; i < fInput ; i++){
381     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
382     if ( !sdigits->GetEntriesFast() )
383       continue ; 
384     Int_t curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
385     if(curNext < nextSig) 
386       nextSig = curNext ;
387   }
388   
389   TArrayI index(fInput) ;
390   index.Reset() ;  //Set all indexes to zero
391   
392   AliPHOSDigit * digit ;
393   AliPHOSDigit * curSDigit ;
394   
395 //  TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
396   
397   //Put Noise contribution
398   Double_t apdNoise = 0. ;
399   if(toMakeNoise)
400      apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ; 
401
402   Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
403   Int_t idigit= 0;
404   for(Int_t imod=0; imod<5; imod++){
405     if(!isPresent[imod])
406       continue ;
407     Int_t firstAbsId=imod*emcpermod+1 ;
408     Int_t lastAbsId =(imod+1)*emcpermod ; 
409     for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
410       Float_t noise = gRandom->Gaus(0.,apdNoise) ; 
411       new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
412       //look if we have to add signal?
413       digit = dynamic_cast<AliPHOSDigit *>(digits->At(idigit)) ;
414       idigit++ ;
415     
416       if(absID==nextSig){
417       //Add SDigits from all inputs 
418 //      ticks->Clear() ;
419 //      Int_t contrib = 0 ;
420
421 //New Timing model is necessary
422 //      Float_t a = digit->GetEnergy() ;
423 //      Float_t b = TMath::Abs( a / fTimeSignalLength) ;
424 //      //Mark the beginning of the signal
425 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);  
426 //      //Mark the end of the signal     
427 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b); 
428
429 // Calculate time as time of the largest digit
430         Float_t time = digit->GetTime() ;
431         Float_t eTime= digit->GetEnergy() ;
432       
433         //loop over inputs
434         for(Int_t i = 0 ; i < fInput ; i++){
435           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
436             curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;     
437             if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
438               curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
439               curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
440             }
441           }
442           else
443             curSDigit = 0 ;
444           //May be several digits will contribute from the same input
445           while(curSDigit && curSDigit->GetId() == absID){         
446             //Shift primary to separate primaries belonging different inputs
447             Int_t primaryoffset ;
448             if(fManager)
449               primaryoffset = fManager->GetMask(i) ; 
450             else
451               primaryoffset = 10000000*i ;
452             curSDigit->ShiftPrimary(primaryoffset) ;
453           
454 //New Timing model is necessary
455 //        a = curSDigit->GetEnergy() ;
456 //        b = a /fTimeSignalLength ;
457 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);  
458 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
459           if(curSDigit->GetEnergy()>eTime){
460             eTime=curSDigit->GetEnergy() ;
461             time=curSDigit->GetTime() ;
462           }
463             *digit += *curSDigit ;  //add energies
464
465             index[i]++ ;
466             if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
467               curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;   
468             else
469               curSDigit = 0 ;
470           }
471         }
472       
473         //calculate and set time
474 //New Timing model is necessary
475 //      Float_t time = FrontEdgeTime(ticks) ;
476         digit->SetTime(time) ;
477       
478         //Find next signal module
479         nextSig = 200000 ;
480         for(Int_t i = 0 ; i < fInput ; i++){
481           sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
482           Int_t curNext = nextSig ;
483           if(sdigits->GetEntriesFast() > index[i] ){
484             curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
485           }
486           if(curNext < nextSig) nextSig = curNext ;
487         }
488       }
489     }
490   }
491
492   //distretize energy if necessary
493   if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
494     Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
495     for(Int_t i = 0 ; i < nEMC ; i++){                                                                                                       
496       digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
497       digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
498     } 
499   }
500   //Apply decalibration if necessary
501   for(Int_t i = 0 ; i < nEMC ; i++){
502     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
503     Decalibrate(digit) ;
504   }
505  
506   
507 //  ticks->Delete() ;
508 //  delete ticks ;
509   
510   //Now CPV digits (different noise and no timing)
511   Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
512   Int_t nEMCtotal=emcpermod*5 ;
513   Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
514   if(isCPVpresent){  //CPV is present in current geometry
515     for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
516       if(!isPresent[imod])
517         continue ;
518       Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
519       Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
520       for(absID = firstAbsId; absID <= lastAbsId; absID++){
521         Float_t noise = gRandom->Gaus(0., cpvNoise) ; 
522         new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
523         idigit++ ;
524         //look if we have to add signal?
525         if(absID==nextSig){
526           digit = dynamic_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
527           //Add SDigits from all inputs
528           for(Int_t i = 0 ; i < fInput ; i++){
529              if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
530                curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;         
531              else
532                curSDigit = 0 ;
533
534              //May be several digits will contribute from the same input
535              while(curSDigit && curSDigit->GetId() == absID){      
536                //Shift primary to separate primaries belonging different inputs
537                Int_t primaryoffset ;
538                if(fManager)
539                  primaryoffset = fManager->GetMask(i) ; 
540                else
541                  primaryoffset = 10000000*i ;
542                curSDigit->ShiftPrimary(primaryoffset) ;
543
544                //add energies
545                *digit += *curSDigit ;  
546                index[i]++ ;
547                if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
548                  curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;      
549                else
550                  curSDigit = 0 ;
551              }
552           }
553
554           //Find next signal module
555           nextSig = 200000 ;
556           for(Int_t i = 0 ; i < fInput ; i++){
557             sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
558             Int_t curNext = nextSig ;
559             if(sdigits->GetEntriesFast() > index[i] )
560               curNext = dynamic_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
561             if(curNext < nextSig) nextSig = curNext ;
562           }
563       
564         }
565       }
566     }
567   }
568
569   delete sdigArray ; //We should not delete its contents 
570   
571   Int_t relId[4];
572
573   //set amplitudes in bad channels to zero
574
575   for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
576     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
577     geom->AbsToRelNumbering(digit->GetId(),relId);
578     if(relId[1] == 0) // Emc
579       if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.); 
580   }
581
582   //remove digits below thresholds
583   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
584   for(Int_t i = 0 ; i < nEMC ; i++){
585     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
586
587     if(digit->GetEnergy() < emcThreshold){
588       digits->RemoveAt(i) ;
589       continue ;
590     }
591
592     digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
593
594     Float_t tres = TimeResolution(digit->GetEnergy()) ; 
595     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
596   }
597
598   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
599   for(Int_t i = nEMC; i < nCPV ; i++){
600     if( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
601       digits->RemoveAt(i) ;
602   } 
603     
604   digits->Compress() ;  
605   Int_t ndigits = digits->GetEntriesFast() ;
606   
607   //Set indexes in list of digits and make true digitization of the energy
608   for (Int_t i = 0 ; i < ndigits ; i++) { 
609     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ; 
610     digit->SetIndexInList(i) ;     
611     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
612       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
613     }
614   }
615
616 }
617 //____________________________________________________________________________
618 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
619   //Apply calibration
620   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
621
622   //Determine rel.position of the cell absolute ID
623   Int_t relId[4];
624   geom->AbsToRelNumbering(absId,relId);
625   Int_t module=relId[0];
626   Int_t row   =relId[2];
627   Int_t column=relId[3];
628   if(relId[1]==0){ //This Is EMC
629     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
630     return amp*calibration ;
631   }
632   return 0 ;
633 }
634 //____________________________________________________________________________
635 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
636 {
637   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
638
639   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
640
641   //Determine rel.position of the cell absolute ID
642   Int_t relId[4];
643   geom->AbsToRelNumbering(digit->GetId(),relId);
644   Int_t module=relId[0];
645   Int_t row   =relId[2];
646   Int_t column=relId[3];
647   if(relId[1]==0){ //This Is EMC
648     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
649     Float_t energy = digit->GetEnergy()/calibration;
650     digit->SetEnergy(energy); //Now digit measures E in ADC counts
651     Float_t time = digit->GetTime() ;
652     time-=fcdb->GetTimeShiftEmc(module,column,row);
653     digit->SetTime(time) ;
654   }
655 }
656 //____________________________________________________________________________
657 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
658   //Apply time calibration
659   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
660
661   //Determine rel.position of the cell absolute ID
662   Int_t relId[4];
663   geom->AbsToRelNumbering(absId,relId);
664   Int_t module=relId[0];
665   Int_t row   =relId[2];
666   Int_t column=relId[3];
667   time += fcdb->GetTimeShiftEmc(module,column,row);
668   return time ;
669 }
670 //____________________________________________________________________________
671 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
672 {
673   // Returns digitized value of the CPV charge in a pad absId
674
675   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
676
677   //Determine rel.position of the cell absId
678   Int_t relId[4];
679   geom->AbsToRelNumbering(absId,relId);
680   Int_t module=relId[0];
681   Int_t row   =relId[2];
682   Int_t column=relId[3];
683   
684   Int_t channel = 0;
685   
686   if(absId > fEmcCrystals){ //digitize CPV only
687
688     //reading calibration data for cell absId.
689     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
690     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
691
692     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;       
693     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
694     if(channel > nMax ) channel = nMax ;
695   }
696   return channel ;
697 }
698
699 //____________________________________________________________________________
700 void AliPHOSDigitizer::Exec(Option_t *option) 
701
702   // Steering method to process digitization for events
703   // in the range from fFirstEvent to fLastEvent.
704   // This range is optionally set by SetEventRange().
705   // if fLastEvent=-1, then process events until the end.
706   // by default fLastEvent = fFirstEvent (process only one event)
707
708   if (!fInit) { // to prevent overwrite existing file
709     AliError(Form("Give a version name different from %s", 
710                   fEventFolderName.Data() )) ;
711     return ;
712   }   
713
714   if (strstr(option,"print")) {
715     Print();
716     return ; 
717   }
718   
719  if(strstr(option,"tim"))
720      gBenchmark->Start("PHOSDigitizer");
721   
722   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
723   AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
724
725   // Post Digitizer to the white board
726   phosLoader->PostDigitizer(this) ;
727   
728   if (fLastEvent == -1) 
729     fLastEvent = rl->GetNumberOfEvents() - 1 ;
730   else if (fManager) 
731     fLastEvent = fFirstEvent ; 
732  
733   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
734   
735   Int_t ievent ;
736
737   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
738     fEventCounter++ ; 
739
740     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
741
742     WriteDigits() ;
743
744     if(strstr(option,"deb"))
745       PrintDigits(option);
746     
747     //increment the total number of Digits per run 
748     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;  
749  }
750  
751   phosLoader->CleanDigitizer();
752
753   if(strstr(option,"tim")){
754     gBenchmark->Stop("PHOSDigitizer");
755     TString message ; 
756     message = "  took %f seconds for Digitizing %f seconds per event\n" ; 
757     AliInfo(Form( message.Data(), 
758          gBenchmark->GetCpuTime("PHOSDigitizer"), 
759          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents )); 
760   } 
761 }
762 //____________________________________________________________________________ 
763 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
764   //calculate TOF resolution using beam-test resutls
765   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
766   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
767   return TMath::Sqrt(a*a+b*b/e) ;
768 }
769
770 ////____________________________________________________________________________ 
771 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
772 //{
773 //  // Returns the shortest time among all time ticks
774 //
775 //  ticks->Sort() ; //Sort in accordance with times of ticks
776 //  TIter it(ticks) ;
777 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
778 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
779 //
780 //  AliPHOSTick * t ;  
781 //  while((t=(AliPHOSTick*) it.Next())){
782 //    if(t->GetTime() < time)  //This tick starts before crossing
783 //      *ctick+=*t ;
784 //    else
785 //      return time ;
786 //
787 //    time = ctick->CrossingTime(fTimeThreshold) ;    
788 //  }
789 //  return time ;
790 //}
791
792 //____________________________________________________________________________ 
793 Bool_t AliPHOSDigitizer::Init()
794 {
795   // Makes all memory allocations
796   fInit = kTRUE ; 
797
798   AliPHOSGeometry *geom;
799   if (!(geom = AliPHOSGeometry::GetInstance())) 
800         geom = AliPHOSGeometry::GetInstance("IHEP","");
801 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
802
803   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
804   
805   fFirstEvent = 0 ; 
806   fLastEvent = fFirstEvent ; 
807   if (fManager) 
808     fInput = fManager->GetNinputs() ; 
809   else 
810     fInput           = 1 ;
811   
812   fInputFileNames  = new TString[fInput] ; 
813   fEventNames      = new TString[fInput] ; 
814   fInputFileNames[0] = GetTitle() ; 
815   fEventNames[0]     = fEventFolderName.Data() ; 
816   Int_t index ; 
817   for (index = 1 ; index < fInput ; index++) {
818     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0); 
819     TString tempo = fManager->GetInputFolderName(index) ;
820     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
821   }
822
823   //to prevent cleaning of this object while GetEvent is called
824   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
825   if(!rl){
826     rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ; 
827   }
828   AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
829   phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
830
831   return fInit ; 
832 }
833
834 //____________________________________________________________________________ 
835 void AliPHOSDigitizer::InitParameters()
836 {
837   // Set initial parameters Digitizer
838
839   fDigitsInRun  = 0 ; 
840   SetEventRange(0,-1) ;
841     
842 }
843
844 //__________________________________________________________________
845 void AliPHOSDigitizer::Print(const Option_t *)const 
846 {
847   // Print Digitizer's parameters
848   AliInfo(Form("\n------------------- %s -------------", GetName() )) ; 
849   if( strcmp(fEventFolderName.Data(), "") != 0 ){
850     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
851     
852     Int_t nStreams ; 
853     if (fManager) 
854       nStreams =  GetNInputStreams() ;
855     else 
856       nStreams = fInput ; 
857     
858     Int_t index = 0 ;  
859     for (index = 0 ; index < nStreams ; index++) {  
860       TString tempo(fEventNames[index]) ; 
861       tempo += index ;
862       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ; 
863     }
864     
865  //   printf("\nWith following parameters:\n") ;
866  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ; 
867  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;  
868  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ; 
869  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;  
870     printf(" ---------------------------------------------------\n") ;   
871   }
872   else
873     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
874   
875 }
876
877 //__________________________________________________________________
878  void AliPHOSDigitizer::PrintDigits(Option_t * option)
879 {
880   // Print a table of digits
881   
882
883   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
884   AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
885   TClonesArray * digits = phosLoader->Digits() ; 
886   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
887   
888   AliInfo(Form("%d", digits->GetEntriesFast())) ; 
889   printf("\nevent %d", gAlice->GetEvNumber()) ;
890   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;  
891
892  
893   if(strstr(option,"all")||strstr(option,"EMC")){  
894     //loop over digits
895     AliPHOSDigit * digit;
896     printf("\nEMC digits (with primaries):\n")  ;
897     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
898     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
899     Int_t index ;
900     for (index = 0 ; (index < digits->GetEntriesFast()) && 
901            (dynamic_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
902       digit = (AliPHOSDigit * )  digits->At(index) ;
903       if(digit->GetNprimary() == 0) 
904         continue;
905 //       printf("%6d  %8d    %6.5e %4d      %2d :",
906 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
907       printf("%6d  %.4f    %6.5e %4d      %2d :",
908               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
909       Int_t iprimary;
910       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
911         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
912       }    
913       printf("\n") ;
914     }
915   }
916   
917   if(strstr(option,"all")||strstr(option,"CPV")){
918     
919     //loop over CPV digits
920     AliPHOSDigit * digit;
921     printf("\nCPV digits (with primaries):\n")  ;
922     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
923     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
924     Int_t index ;
925     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
926       digit = (AliPHOSDigit * )  digits->At(index) ;
927       if(digit->GetId() > maxEmc){
928         printf("%6d  %8d    %4d      %2d :",
929                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;  
930         Int_t iprimary;
931         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
932           printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
933         }    
934         printf("\n") ;
935       }
936     }
937   }
938
939 }
940
941 //__________________________________________________________________
942 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
943 {  // Calculates the time signal generated by noise
944   //PH  Info("TimeOfNoise", "Change me") ; 
945   return gRandom->Rndm() * 1.28E-5;
946 }
947
948 //__________________________________________________________________
949 void AliPHOSDigitizer::Unload() 
950 {  
951   
952   Int_t i ; 
953   for(i = 1 ; i < fInput ; i++){
954     TString tempo(fEventNames[i]) ; 
955     tempo += i ;
956     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ; 
957     AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
958     phosLoader->UnloadSDigits() ; 
959   }
960   
961   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
962   AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
963   phosLoader->UnloadDigits() ; 
964 }
965
966 //____________________________________________________________________________
967 void AliPHOSDigitizer::WriteDigits()
968 {
969
970   // Makes TreeD in the output file. 
971   // Check if branch already exists: 
972   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
973   //      already existing branches. 
974   //   else creates branch with Digits, named "PHOS", title "...",
975   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
976   //      and names of files, from which digits are made.
977
978   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
979   AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
980  
981   const TClonesArray * digits = phosLoader->Digits() ; 
982   TTree * treeD = phosLoader->TreeD();
983   if(!treeD){
984     phosLoader->MakeTree("D");
985     treeD = phosLoader->TreeD();
986   }
987   
988   // -- create Digits branch
989   Int_t bufferSize = 32000 ;    
990   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
991   digitsBranch->SetTitle(fEventFolderName);
992   digitsBranch->Fill() ;
993   
994   phosLoader->WriteDigits("OVERWRITE");
995   phosLoader->WriteDigitizer("OVERWRITE");
996
997   Unload() ; 
998
999 }