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