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