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