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