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