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