cc874cc923199456538c57592368912bff6af8e7
[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$
21  * Revision 1.98  2007/09/26 14:22:17  cvetan
22  * 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.
23  *
24  * Revision 1.97  2007/08/07 14:12:03  kharlov
25  * Quality assurance added (Yves Schutz)
26  *
27  * Revision 1.96  2007/04/28 10:43:36  policheh
28  * Dead channels simulation: digit energy sets to 0.
29  *
30  * Revision 1.95  2007/04/10 07:20:52  kharlov
31  * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
32  *
33  * Revision 1.94  2007/02/01 10:34:47  hristov
34  * Removing warnings on Solaris x86
35  *
36  * Revision 1.93  2006/10/17 13:17:01  kharlov
37  * Replace AliInfo by AliDebug
38  *
39  * Revision 1.92  2006/08/28 10:01:56  kharlov
40  * Effective C++ warnings fixed (Timur Pocheptsov)
41  *
42  * Revision 1.91  2006/04/29 20:25:30  hristov
43  * Decalibration is implemented (Yu.Kharlov)
44  *
45  * Revision 1.90  2006/04/22 10:30:17  hristov
46  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
47  *
48  * Revision 1.89  2006/04/11 15:22:59  hristov
49  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
50  *
51  * Revision 1.88  2006/03/13 14:05:43  kharlov
52  * Calibration objects for EMC and CPV
53  *
54  * Revision 1.87  2005/08/24 15:33:49  kharlov
55  * Calibration data for raw digits
56  *
57  * Revision 1.86  2005/07/12 20:07:35  hristov
58  * Changes needed to run simulation and reconstrruction in the same AliRoot session
59  *
60  * Revision 1.85  2005/05/28 14:19:04  schutz
61  * Compilation warnings fixed by T.P.
62  *
63  */
64
65 //_________________________________________________________________________
66 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute) 
67 //////////////////////////////////////////////////////////////////////////////
68 // This TTask performs digitization of Summable digits (in the PHOS case it is just
69 // the sum of contributions from all primary particles into a given cell). 
70 // In addition it performs mixing of summable digits from different events.
71 // The name of the TTask is also the title of the branch that will contain 
72 // the created SDigits
73 // The title of the TTAsk is the name of the file that contains the hits from
74 // which the SDigits are created
75 //
76 // For each event two branches are created in TreeD:
77 //   "PHOS" - list of digits
78 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
79 //
80 // Note, that one can set a title for new digits branch, and repeat digitization with
81 // another set of parameters.
82 //
83 // Use case:
84 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
85 // root[1] d->ExecuteTask()             
86 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
87 //                       //Digitizes SDigitis in all events found in file galice.root 
88 //
89 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;  
90 //                       // Will read sdigits from galice1.root
91 // root[3] d1->MixWith("galice2.root")       
92 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
93 //                       // Reads another set of sdigits from galice2.root
94 // root[3] d1->MixWith("galice3.root")       
95 //                       // Reads another set of sdigits from galice3.root
96 // root[4] d->ExecuteTask("deb timing")    
97 //                       // Reads SDigits from files galice1.root, galice2.root ....
98 //                       // mixes them and stores produced Digits in file galice1.root          
99 //                       // deb - prints number of produced digits
100 //                       // deb all - prints list of produced digits
101 //                       // timing  - prints time used for digitization
102 //
103
104 // --- ROOT system ---
105 #include "TTree.h"
106 #include "TSystem.h"
107 #include "TBenchmark.h"
108 #include "TRandom.h"
109
110 // --- Standard library ---
111
112 // --- AliRoot header files ---
113 #include "AliLog.h"
114 #include "AliRunDigitizer.h"
115 #include "AliPHOSDigit.h"
116 #include "AliPHOSGetter.h"
117 #include "AliPHOSDigitizer.h"
118 #include "AliPHOSSDigitizer.h"
119 #include "AliPHOSGeometry.h"
120 #include "AliPHOSTick.h"
121 #include "AliPHOSQualAssDataMaker.h" 
122
123 ClassImp(AliPHOSDigitizer)
124
125
126 //____________________________________________________________________________ 
127 AliPHOSDigitizer::AliPHOSDigitizer() :
128   AliDigitizer("",""),
129   fDefaultInit(kTRUE),
130   fDigitsInRun(0),
131   fInit(kFALSE),
132   fInput(0),
133   fInputFileNames(0x0),
134   fEventNames(0x0),
135   fEmcCrystals(0),
136   fPinNoise(0.f),
137   fEMCDigitThreshold(0.f),
138   fCPVNoise(0.f),
139   fCPVDigitThreshold(0.f),
140   fTimeResolution(0.f),
141   fTimeThreshold(0.f),
142   fTimeSignalLength(0.f),
143   fADCchanelEmc(0.f),
144   fADCpedestalEmc(0.f),
145   fNADCemc(0),
146   fADCchanelCpv(0.f),
147   fADCpedestalCpv(0.f),
148   fNADCcpv(0),
149   fEventFolderName(""),
150   fFirstEvent(0),
151   fLastEvent(0), 
152   fQADM (0x0), 
153   fEventCounter(0)
154 {
155   // ctor
156   InitParameters() ; 
157   fManager = 0 ;                     // We work in the standalong mode
158 }
159
160 //____________________________________________________________________________ 
161 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName, 
162                                    TString eventFolderName):
163   AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName), 
164   fDefaultInit(kFALSE),
165   fDigitsInRun(0),
166   fInit(kFALSE),
167   fInput(0),
168   fInputFileNames(0x0),
169   fEventNames(0x0),
170   fEmcCrystals(0),
171   fPinNoise(0.f),
172   fEMCDigitThreshold(0.f),
173   fCPVNoise(0.f),
174   fCPVDigitThreshold(0.f),
175   fTimeResolution(0.f),
176   fTimeThreshold(0.f),
177   fTimeSignalLength(0.f),
178   fADCchanelEmc(0.f),
179   fADCpedestalEmc(0.f),
180   fNADCemc(0),
181   fADCchanelCpv(0.f),
182   fADCpedestalCpv(0.f),
183   fNADCcpv(0),
184   fEventFolderName(eventFolderName),
185   fFirstEvent(0),
186   fLastEvent(0), 
187   fQADM (0x0), 
188   fEventCounter(0)
189 {
190   // ctor
191   InitParameters() ; 
192   Init() ;
193   fDefaultInit = kFALSE ; 
194   fManager = 0 ;                     // We work in the standalong mode
195   //Initialize the quality assurance data maker only once
196   fQADM = new AliPHOSQualAssDataMaker() ;  
197   //FIXME: get the run number
198   Int_t run = 0 ;
199   //EMXIF       
200    GetQualAssDataMaker()->Init(AliQualAss::kDIGITS, run, fgkCycles) ;    
201    GetQualAssDataMaker()->StartOfCycle(AliQualAss::kDIGITS) ;    
202 }
203
204 //____________________________________________________________________________ 
205 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) : 
206   AliDigitizer(d),
207   fDefaultInit(d.fDefaultInit),
208   fDigitsInRun(d.fDigitsInRun),
209   fInit(d.fInit),
210   fInput(d.fInput),
211   fInputFileNames(0x0),//?
212   fEventNames(0x0),//?
213   fEmcCrystals(d.fEmcCrystals),
214   fPinNoise(d.fPinNoise),
215   fEMCDigitThreshold(d.fEMCDigitThreshold),
216   fCPVNoise(d.fCPVNoise),
217   fCPVDigitThreshold(d.fCPVDigitThreshold),
218   fTimeResolution(d.fTimeResolution),
219   fTimeThreshold(d.fTimeThreshold),
220   fTimeSignalLength(d.fTimeSignalLength),
221   fADCchanelEmc(d.fADCchanelEmc),
222   fADCpedestalEmc(d.fADCpedestalEmc),
223   fNADCemc(d.fNADCemc),
224   fADCchanelCpv(d.fADCchanelCpv),
225   fADCpedestalCpv(d.fADCpedestalCpv),
226   fNADCcpv(d.fNADCcpv),
227   fEventFolderName(d.fEventFolderName),
228   fFirstEvent(d.fFirstEvent),
229   fLastEvent(d.fLastEvent), 
230   fQADM (d.fQADM), 
231   fEventCounter(0)
232
233 {
234   // copyy ctor 
235   SetName(d.GetName()) ; 
236   SetTitle(d.GetTitle()) ; 
237 //Initialize the quality assurance data maker only once
238  //FIXME: get the run number
239   Int_t run = 0 ;
240   //EMXIF       
241   GetQualAssDataMaker()->Init(AliQualAss::kDIGITS, run, fgkCycles) ;    
242   GetQualAssDataMaker()->StartOfCycle(AliQualAss::kDIGITS) ;    
243 }
244
245 //____________________________________________________________________________ 
246 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
247   AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
248   fDefaultInit(kFALSE),
249   fDigitsInRun(0),
250   fInit(kFALSE),
251   fInput(0),
252   fInputFileNames(0x0),
253   fEventNames(0x0),
254   fEmcCrystals(0),
255   fPinNoise(0.f),
256   fEMCDigitThreshold(0.f),
257   fCPVNoise(0.f),
258   fCPVDigitThreshold(0.f),
259   fTimeResolution(0.f),
260   fTimeThreshold(0.f),
261   fTimeSignalLength(0.f),
262   fADCchanelEmc(0.f),
263   fADCpedestalEmc(0.f),
264   fNADCemc(0),
265   fADCchanelCpv(0.f),
266   fADCpedestalCpv(0.f),
267   fNADCcpv(0),
268   fEventFolderName(fManager->GetInputFolderName(0)),
269   fFirstEvent(0),
270   fLastEvent(0), 
271   fQADM (0x0), 
272   fEventCounter(0)
273
274 {
275   // ctor Init() is called by RunDigitizer
276   fManager = rd ; 
277   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
278   InitParameters() ; 
279   fDefaultInit = kFALSE ; 
280 //Initialize the quality assurance data maker only once
281   fQADM = new AliPHOSQualAssDataMaker() ;  
282  //FIXME: get the run number
283   Int_t run = 0 ;
284   //EMXIF       
285   GetQualAssDataMaker()->Init(AliQualAss::kDIGITS, run) ;    
286   GetQualAssDataMaker()->StartOfCycle(AliQualAss::kDIGITS) ;    
287 }
288
289 //____________________________________________________________________________ 
290   AliPHOSDigitizer::~AliPHOSDigitizer()
291 {
292   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
293
294   // Clean Digitizer from the white board
295   gime->PhosLoader()->CleanDigitizer() ;
296   // dtor
297   delete [] fInputFileNames ; 
298   delete [] fEventNames ; 
299
300   delete fQADM ; 
301
302 }
303
304 //____________________________________________________________________________
305 void AliPHOSDigitizer::Digitize(Int_t event) 
306
307   
308   // Makes the digitization of the collected summable digits.
309   //  It first creates the array of all PHOS modules
310   //  filled with noise (different for EMC, and CPV) and
311   //  then adds contributions from SDigits. 
312   // This design avoids scanning over the list of digits to add 
313   // contribution to new SDigits only.
314
315   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ; 
316   Int_t ReadEvent = event ; 
317   if (fManager) 
318     ReadEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ; 
319   AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
320                   ReadEvent, GetTitle(), fEventFolderName.Data())) ; 
321   gime->Event(ReadEvent, "S") ;
322   TClonesArray * digits = gime->Digits() ; 
323   digits->Clear() ;
324
325   const AliPHOSGeometry *geom = gime->PHOSGeometry() ; 
326   //Making digits with noise, first EMC
327   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
328   
329   Int_t nCPV ;
330   Int_t absID ;
331   
332   nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * geom->GetNModules() ;
333   
334   digits->Expand(nCPV) ;
335
336   // get first the sdigitizer from the tasks list 
337   if ( !gime->SDigitizer() ) 
338     gime->LoadSDigitizer();
339   AliPHOSSDigitizer * sDigitizer = gime->SDigitizer(); 
340
341   if ( !sDigitizer )
342     AliFatal(Form("SDigitizer with name %s %s not found", 
343                   GetTitle(), fEventFolderName.Data() )) ; 
344
345   //take all the inputs to add together and load the SDigits
346   TObjArray * sdigArray = new TObjArray(fInput) ;
347   sdigArray->AddAt(gime->SDigits(), 0) ;
348   Int_t i ;
349   for(i = 1 ; i < fInput ; i++){
350     TString tempo(fEventNames[i]) ; 
351     tempo += i ;
352     AliPHOSGetter * gime1 = AliPHOSGetter::Instance(fInputFileNames[i], tempo) ;
353     if (fManager) 
354       ReadEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ; 
355     AliInfo(Form("Adding event %d from input stream %d %s %s", 
356                  ReadEvent, i, fInputFileNames[i].Data(), tempo.Data())) ; 
357     gime1->Event(ReadEvent,"S");
358     sdigArray->AddAt(gime1->SDigits(), i) ;
359   }
360
361   //Find the first crystal with signal
362   Int_t nextSig = 200000 ; 
363   TClonesArray * sdigits ;  
364   for(i = 0 ; i < fInput ; i++){
365     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
366     if ( !sdigits->GetEntriesFast() )
367       continue ; 
368     Int_t curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
369     if(curNext < nextSig) 
370       nextSig = curNext ;
371   }
372   
373   TArrayI index(fInput) ;
374   index.Reset() ;  //Set all indexes to zero
375   
376   AliPHOSDigit * digit ;
377   AliPHOSDigit * curSDigit ;
378   
379   TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
380   
381   //Put Noise contribution
382   for(absID = 1 ; absID <= nEMC ; absID++){
383     Float_t noise = gRandom->Gaus(0., fPinNoise) ; 
384     // YVK: do not digitize amplitudes for EMC
385 //     new((*digits)[absID-1]) AliPHOSDigit( -1, absID, sDigitizer->Digitize(noise), TimeOfNoise() ) ;
386     new((*digits)[absID-1]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
387     //look if we have to add signal?
388     digit = dynamic_cast<AliPHOSDigit *>(digits->At(absID-1)) ;
389     
390     if(absID==nextSig){
391       //Add SDigits from all inputs 
392       ticks->Clear() ;
393       Int_t contrib = 0 ;
394       Float_t a = digit->GetEnergy() ;
395       Float_t b = TMath::Abs( a / fTimeSignalLength) ;
396       //Mark the beginning of the signal
397       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);  
398       //Mark the end of the signal     
399       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b); 
400       
401       //loop over inputs
402       for(i = 0 ; i < fInput ; i++){
403         if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
404           curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;       
405         else
406           curSDigit = 0 ;
407         //May be several digits will contribute from the same input
408         while(curSDigit && curSDigit->GetId() == absID){           
409           //Shift primary to separate primaries belonging different inputs
410           Int_t primaryoffset ;
411           if(fManager)
412             primaryoffset = fManager->GetMask(i) ; 
413           else
414             primaryoffset = 10000000*i ;
415           curSDigit->ShiftPrimary(primaryoffset) ;
416           
417           a = curSDigit->GetEnergy() ;
418           b = a /fTimeSignalLength ;
419           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);  
420           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
421           
422           *digit += *curSDigit ;  //add energies
423
424           index[i]++ ;
425           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
426             curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;     
427           else
428             curSDigit = 0 ;
429         }
430       }
431       
432       //calculate and set time
433       Float_t time = FrontEdgeTime(ticks) ;
434       digit->SetTime(time) ;
435       
436       //Find next signal module
437       nextSig = 200000 ;
438       for(i = 0 ; i < fInput ; i++){
439         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
440         Int_t curNext = nextSig ;
441         if(sdigits->GetEntriesFast() > index[i] ){
442           curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
443         }
444         if(curNext < nextSig) nextSig = curNext ;
445       }
446     }
447   }
448   
449   ticks->Delete() ;
450   delete ticks ;
451   
452   //Now CPV digits (different noise and no timing)
453   for(absID = nEMC+1; absID <= nCPV; absID++){
454     Float_t noise = gRandom->Gaus(0., fCPVNoise) ; 
455     new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
456     //look if we have to add signal?
457     if(absID==nextSig){
458       digit = dynamic_cast<AliPHOSDigit *>(digits->At(absID-1)) ;
459       //Add SDigits from all inputs
460       for(i = 0 ; i < fInput ; i++){
461         if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
462           curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;      
463         else
464           curSDigit = 0 ;
465
466         //May be several digits will contribute from the same input
467         while(curSDigit && curSDigit->GetId() == absID){           
468           //Shift primary to separate primaries belonging different inputs
469           Int_t primaryoffset ;
470           if(fManager)
471             primaryoffset = fManager->GetMask(i) ; 
472           else
473             primaryoffset = 10000000*i ;
474           curSDigit->ShiftPrimary(primaryoffset) ;
475
476           //add energies
477           *digit += *curSDigit ;  
478           index[i]++ ;
479           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
480             curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;   
481           else
482             curSDigit = 0 ;
483         }
484       }
485
486       //Find next signal module
487       nextSig = 200000 ;
488       for(i = 0 ; i < fInput ; i++){
489         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
490         Int_t curNext = nextSig ;
491         if(sdigits->GetEntriesFast() > index[i] )
492           curNext = dynamic_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
493         if(curNext < nextSig) nextSig = curNext ;
494       }
495       
496     }
497   }
498
499   delete sdigArray ; //We should not delete its contents 
500   
501   //remove digits below thresholds
502   for(i = 0 ; i < nEMC ; i++){
503     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
504     DecalibrateEMC(digit);
505     if(digit->GetEnergy() < fEMCDigitThreshold)
506       digits->RemoveAt(i) ;
507     else
508       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
509   }
510
511
512   for(i = nEMC; i < nCPV ; i++)
513 //     if( sDigitizer->Calibrate( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetAmp() ) < fCPVDigitThreshold )
514     if( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < fCPVDigitThreshold )
515       digits->RemoveAt(i) ;
516     
517   digits->Compress() ;  
518   
519   Int_t ndigits = digits->GetEntriesFast() ;
520   digits->Expand(ndigits) ;
521
522   //Set indexes in list of digits and make true digitization of the energy
523   for (i = 0 ; i < ndigits ; i++) { 
524     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ; 
525     digit->SetIndexInList(i) ;     
526     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
527       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
528     }
529   }
530
531   Int_t relId[4];
532
533   //set amplitudes in bad channels to zero
534   for(i = 0 ; i <digits->GetEntries(); i++){
535     digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
536     gime->PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
537     if(relId[1] == 0) // Emc
538       if(gime->CalibData()->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.); 
539   }
540
541 }
542
543 //____________________________________________________________________________
544 void AliPHOSDigitizer::DecalibrateEMC(AliPHOSDigit *digit)
545 {
546   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
547
548   AliPHOSGetter* gime = AliPHOSGetter::Instance();
549
550   if(!gime->CalibData()) {
551     AliPHOSCalibData* cdb = new AliPHOSCalibData(-1);
552     gime->SetCalibData(cdb);
553   }
554
555   //Determine rel.position of the cell absolute ID
556   Int_t relId[4];
557   gime->PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
558   Int_t module=relId[0];
559   Int_t row   =relId[2];
560   Int_t column=relId[3];
561   Float_t decalibration = gime->CalibData()->GetADCchannelEmc(module,column,row);
562   Float_t energy = digit->GetEnergy() / decalibration;
563   digit->SetEnergy(energy);
564 }
565 //____________________________________________________________________________
566 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
567 {
568   // Returns digitized value of the CPV charge in a pad absId
569
570   AliPHOSGetter* gime = AliPHOSGetter::Instance();
571
572   if(!gime->CalibData()) {
573     AliPHOSCalibData* cdb = new AliPHOSCalibData(-1); // use AliCDBManager's run number
574     gime->SetCalibData(cdb);
575   }
576
577   //Determine rel.position of the cell absId
578   Int_t relId[4];
579   gime->PHOSGeometry()->AbsToRelNumbering(absId,relId);
580   Int_t module=relId[0];
581   Int_t row   =relId[2];
582   Int_t column=relId[3];
583   
584   Int_t channel = 0;
585   
586   if(absId > fEmcCrystals){ //digitize CPV only
587
588     //reading calibration data for cell absId.
589     //If no calibration DB found, accept default values.
590
591     if(gime->CalibData()) {
592       fADCpedestalCpv = gime->CalibData()->GetADCpedestalCpv(module,column,row);
593       fADCchanelCpv   = gime->CalibData()->GetADCchannelCpv( module,column,row);
594     }
595
596     channel = (Int_t) TMath::Ceil((charge - fADCpedestalCpv)/fADCchanelCpv) ;       
597     if(channel > fNADCcpv ) channel =  fNADCcpv ;
598   }
599   return channel ;
600 }
601
602 //____________________________________________________________________________
603 void AliPHOSDigitizer::Exec(Option_t *option) 
604
605   // Steering method to process digitization for events
606   // in the range from fFirstEvent to fLastEvent.
607   // This range is optionally set by SetEventRange().
608   // if fLastEvent=-1, then process events until the end.
609   // by default fLastEvent = fFirstEvent (process only one event)
610
611   if (!fInit) { // to prevent overwrite existing file
612     AliError(Form("Give a version name different from %s", 
613                   fEventFolderName.Data() )) ;
614     return ;
615   }   
616
617   if (strstr(option,"print")) {
618     Print();
619     return ; 
620   }
621   
622   if(strstr(option,"tim"))
623     gBenchmark->Start("PHOSDigitizer");
624   
625   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
626
627   // Post Digitizer to the white board
628   gime->PostDigitizer(this) ;
629   
630   if (fLastEvent == -1) 
631     fLastEvent = gime->MaxEvent() - 1 ;
632   else if (fManager) 
633     fLastEvent = fFirstEvent ; 
634  
635   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
636   
637   Int_t ievent ;
638
639   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
640     fEventCounter++ ; 
641     gime->Event(ievent,"S") ;
642
643     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
644
645     //makes the quality assurance data
646     if (GetQualAssDataMaker()->IsCycleDone() ) {
647       GetQualAssDataMaker()->EndOfCycle(AliQualAss::kDIGITS) ; 
648           GetQualAssDataMaker()->StartOfCycle(AliQualAss::kDIGITS) ; 
649     }
650     GetQualAssDataMaker()->Exec(AliQualAss::kDIGITS, gime->Digits()) ; 
651         GetQualAssDataMaker()->Increment() ;
652                                    
653     WriteDigits() ;
654
655     if(strstr(option,"deb"))
656       PrintDigits(option);
657     
658     //increment the total number of Digits per run 
659     fDigitsInRun += gime->Digits()->GetEntriesFast() ;  
660  }
661  
662    //Write the quality assurance data only after the last event 
663  if ( fEventCounter == gime->MaxEvent() ) { 
664     GetQualAssDataMaker()->EndOfCycle(AliQualAss::kDIGITS) ;    
665         GetQualAssDataMaker()->Finish(AliQualAss::kDIGITS) ;
666   }
667   
668   gime->PhosLoader()->CleanDigitizer();
669
670   if(strstr(option,"tim")){
671     gBenchmark->Stop("PHOSDigitizer");
672     TString message ; 
673     message = "  took %f seconds for Digitizing %f seconds per event\n" ; 
674     AliInfo(Form( message.Data(), 
675          gBenchmark->GetCpuTime("PHOSDigitizer"), 
676          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents )); 
677   } 
678 }
679
680 //____________________________________________________________________________ 
681 Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
682 {
683   // Returns the shortest time among all time ticks
684
685   ticks->Sort() ; //Sort in accordance with times of ticks
686   TIter it(ticks) ;
687   AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
688   Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
689
690   AliPHOSTick * t ;  
691   while((t=(AliPHOSTick*) it.Next())){
692     if(t->GetTime() < time)  //This tick starts before crossing
693       *ctick+=*t ;
694     else
695       return time ;
696
697     time = ctick->CrossingTime(fTimeThreshold) ;    
698   }
699   return time ;
700 }
701
702 //____________________________________________________________________________ 
703 Bool_t AliPHOSDigitizer::Init()
704 {
705   // Makes all memory allocations
706   fInit = kTRUE ; 
707   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ; 
708   if ( gime == 0 ) {
709     AliFatal(Form("Could not obtain the Getter object for file %s and event %s !", 
710                   GetTitle(), fEventFolderName.Data()));  
711     return kFALSE;
712   } 
713   
714   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
715
716   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
717   
718   TString opt("Digits") ; 
719   if(gime->VersionExists(opt) ) { 
720     AliError(Form("Give a version name different from %s", 
721                   fEventFolderName.Data() )) ;
722     fInit = kFALSE ; 
723   }
724
725   fFirstEvent = 0 ; 
726   fLastEvent = fFirstEvent ; 
727   if (fManager) 
728     fInput = fManager->GetNinputs() ; 
729   else 
730     fInput           = 1 ;
731   
732   fInputFileNames  = new TString[fInput] ; 
733   fEventNames      = new TString[fInput] ; 
734   fInputFileNames[0] = GetTitle() ; 
735   fEventNames[0]     = fEventFolderName.Data() ; 
736   Int_t index ; 
737   for (index = 1 ; index < fInput ; index++) {
738     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0); 
739     TString tempo = fManager->GetInputFolderName(index) ;
740     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
741   }
742
743   //to prevent cleaning of this object while GetEvent is called
744   gime->PhosLoader()->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
745
746   return fInit ; 
747 }
748
749 //____________________________________________________________________________ 
750 void AliPHOSDigitizer::InitParameters()
751 {
752   // Set initial parameters Digitizer
753
754   fPinNoise           = 0.004 ;  // [GeV]
755   fEMCDigitThreshold  = 0.012 ;  // [GeV]
756   fCPVNoise           = 0.01;    // [aux units]
757   fCPVDigitThreshold  = 0.09 ;   // [aux units]
758   fTimeResolution     = 0.5e-9 ; // [sec]
759   fTimeSignalLength   = 1.0e-9 ; // [sec]
760   fDigitsInRun  = 0 ; 
761   fADCchanelEmc = 1.0;        // Coefficient between real and measured energies in EMC
762   fADCpedestalEmc = 0. ;      //
763   fNADCemc = (Int_t) TMath::Power(2,16) ;  // number of channels in EMC ADC
764
765   fADCchanelCpv = 0.0012 ;          // width of one ADC channel in CPV 'popugais'
766   fADCpedestalCpv = 0.012 ;         // 
767   fNADCcpv = (Int_t) TMath::Power(2,12);      // number of channels in CPV ADC
768
769 //   fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
770   fTimeThreshold = 0.001 ; // [GeV]
771   SetEventRange(0,-1) ;
772     
773 }
774
775 //__________________________________________________________________
776 void AliPHOSDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
777 {
778   // Allows to produce digits by superimposing background and signal event.
779   // It is assumed, that headers file with SIGNAL events is opened in 
780   // the constructor. 
781   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed 
782   // Thus we avoid writing (changing) huge and expensive 
783   // backgound files: all output will be writen into SIGNAL, i.e. 
784   // opened in constructor file. 
785   //
786   // One can open as many files to mix with as one needs.
787   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
788   // can be mixed.
789
790   if( strcmp(fEventFolderName, "") == 0 )
791     Init() ;
792
793   if(fManager){
794     Warning("MixWith", "Cannot use this method with AliRunDigitizer\n" ) ;
795     return ;
796   }
797   // looking for file which contains AliRun
798   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
799     AliError(Form("File %s does not exist!", alirunFileName.Data())) ;
800     return ; 
801   }
802   // looking for the file which contains SDigits
803   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
804   TString fileName( gime->GetSDigitsFileName() ) ; 
805     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
806       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
807     if ( (gSystem->AccessPathName(fileName)) ) { 
808       AliError(Form("The file %s does not exist!", fileName.Data())) ;
809       return ;
810     }
811     // need to increase the arrays
812     TString tempo = fInputFileNames[fInput-1] ; 
813     delete [] fInputFileNames ; 
814     fInputFileNames = new TString[fInput+1] ; 
815     fInputFileNames[fInput-1] = tempo ; 
816  
817     tempo = fEventNames[fInput-1] ; 
818     delete [] fEventNames ; 
819     fEventNames = new TString[fInput+1] ; 
820     fEventNames[fInput-1] = tempo ; 
821
822     fInputFileNames[fInput] = alirunFileName ; 
823     fEventNames[fInput]     = eventFolderName ;
824     fInput++ ;
825 }
826
827 //__________________________________________________________________
828 void AliPHOSDigitizer::Print(const Option_t *)const 
829 {
830   // Print Digitizer's parameters
831   AliInfo(Form("\n------------------- %s -------------", GetName() )) ; 
832   if( strcmp(fEventFolderName.Data(), "") != 0 ){
833     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
834     
835     Int_t nStreams ; 
836     if (fManager) 
837       nStreams =  GetNInputStreams() ;
838     else 
839       nStreams = fInput ; 
840     
841     Int_t index = 0 ;  
842     for (index = 0 ; index < nStreams ; index++) {  
843       TString tempo(fEventNames[index]) ; 
844       tempo += index ;
845       AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[index], tempo) ; 
846       TString fileName( gime->GetSDigitsFileName() ) ; 
847       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
848         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
849       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ; 
850     }
851     AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ; 
852     printf("\nWriting digits to %s", gime->GetDigitsFileName().Data()) ;
853     
854     printf("\nWith following parameters:\n") ;
855     printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ; 
856     printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;  
857     printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ; 
858     printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;  
859     printf(" ---------------------------------------------------\n") ;   
860   }
861   else
862     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
863   
864 }
865
866 //__________________________________________________________________
867  void AliPHOSDigitizer::PrintDigits(Option_t * option)
868 {
869   // Print a table of digits
870   
871   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ; 
872   TClonesArray * digits = gime->Digits() ; 
873   
874   AliInfo(Form("%d", digits->GetEntriesFast())) ; 
875   printf("\nevent %d", gAlice->GetEvNumber()) ;
876   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;  
877
878  
879   if(strstr(option,"all")||strstr(option,"EMC")){  
880     //loop over digits
881     AliPHOSDigit * digit;
882     printf("\nEMC digits (with primaries):\n")  ;
883     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
884     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
885     Int_t index ;
886     for (index = 0 ; (index < digits->GetEntriesFast()) && 
887            (dynamic_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
888       digit = (AliPHOSDigit * )  digits->At(index) ;
889       if(digit->GetNprimary() == 0) 
890         continue;
891 //       printf("%6d  %8d    %6.5e %4d      %2d :",
892 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
893       printf("%6d  %.4f    %6.5e %4d      %2d :",
894               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
895       Int_t iprimary;
896       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
897         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
898       }    
899       printf("\n") ;
900     }
901   }
902   
903   if(strstr(option,"all")||strstr(option,"CPV")){
904     
905     //loop over CPV digits
906     AliPHOSDigit * digit;
907     printf("\nCPV digits (with primaries):\n")  ;
908     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
909     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
910     Int_t index ;
911     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
912       digit = (AliPHOSDigit * )  digits->At(index) ;
913       if(digit->GetId() > maxEmc){
914         printf("%6d  %8d    %4d      %2d :",
915                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;  
916         Int_t iprimary;
917         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
918           printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
919         }    
920         printf("\n") ;
921       }
922     }
923   }
924
925 }
926
927 //__________________________________________________________________
928 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
929 {  // Calculates the time signal generated by noise
930   //PH  Info("TimeOfNoise", "Change me") ; 
931   return gRandom->Rndm() * 1.28E-5;
932 }
933
934 //__________________________________________________________________
935 void AliPHOSDigitizer::Unload() 
936 {  
937   
938   Int_t i ; 
939   for(i = 1 ; i < fInput ; i++){
940     TString tempo(fEventNames[i]) ; 
941     tempo += i ;
942     AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileNames[i], tempo) ; 
943     gime->PhosLoader()->UnloadSDigits() ; 
944   }
945   
946   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName) ; 
947   gime->PhosLoader()->UnloadDigits() ; 
948 }
949
950 //____________________________________________________________________________
951 void AliPHOSDigitizer::WriteDigits()
952 {
953
954   // Makes TreeD in the output file. 
955   // Check if branch already exists: 
956   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
957   //      already existing branches. 
958   //   else creates branch with Digits, named "PHOS", title "...",
959   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
960   //      and names of files, from which digits are made.
961
962   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ; 
963   const TClonesArray * digits = gime->Digits() ; 
964   TTree * treeD = gime->TreeD();
965
966   // -- create Digits branch
967   Int_t bufferSize = 32000 ;    
968   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
969   digitsBranch->SetTitle(fEventFolderName);
970   digitsBranch->Fill() ;
971   
972   gime->WriteDigits("OVERWRITE");
973   gime->WriteDigitizer("OVERWRITE");
974
975   Unload() ; 
976
977 }