Fixed memory leaks for #86360: High memory consumption in 2.76TeV p+p RAW reco jobs
[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 TTask 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 TTask is also the title of the branch that will contain 
90 // the created SDigits
91 // The title of the TTAsk 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->ExecuteTask()             
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->ExecuteTask("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 "AliRunDigitizer.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   fManager = 0 ;                     // We work in the standalong mode
170 }
171
172 //____________________________________________________________________________ 
173 AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName, 
174                                    TString eventFolderName):
175   AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), 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   fManager = 0 ;                     // We work in the standalone mode
197   fcdb = new AliPHOSCalibData(-1);
198 }
199
200 //____________________________________________________________________________ 
201 AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) : 
202   AliDigitizer(d),
203   fDefaultInit(d.fDefaultInit),
204   fDigitsInRun(d.fDigitsInRun),
205   fInit(d.fInit),
206   fInput(d.fInput),
207   fInputFileNames(0x0),//?
208   fEventNames(0x0),//?
209   fEmcCrystals(d.fEmcCrystals),
210   fEventFolderName(d.fEventFolderName),
211   fFirstEvent(d.fFirstEvent),
212   fLastEvent(d.fLastEvent), 
213   fcdb (0x0), 
214   fEventCounter(0),
215   fPulse(0),
216   fADCValuesLG(0),
217   fADCValuesHG(0)
218 {
219   // copyy ctor 
220   SetName(d.GetName()) ; 
221   SetTitle(d.GetTitle()) ; 
222   fcdb = new AliPHOSCalibData(-1);
223 }
224
225 //____________________________________________________________________________ 
226 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
227   AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
228   fDefaultInit(kFALSE),
229   fDigitsInRun(0),
230   fInit(kFALSE),
231   fInput(0),
232   fInputFileNames(0x0),
233   fEventNames(0x0),
234   fEmcCrystals(0),
235   fEventFolderName(fManager->GetInputFolderName(0)),
236   fFirstEvent(0),
237   fLastEvent(0), 
238   fcdb (0x0), 
239   fEventCounter(0),
240   fPulse(0),
241   fADCValuesLG(0),
242   fADCValuesHG(0)
243
244 {
245   // ctor Init() is called by RunDigitizer
246   fManager = rd ; 
247   SetTitle(static_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
248   InitParameters() ; 
249   fDefaultInit = kFALSE ; 
250   fcdb = new AliPHOSCalibData(-1);
251 }
252
253 //____________________________________________________________________________ 
254   AliPHOSDigitizer::~AliPHOSDigitizer()
255 {
256   // dtor
257   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
258   if(rl){                                                                                                                               
259     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));    
260     
261     if(phosLoader)
262       phosLoader->CleanDigitizer() ;
263   }
264   
265   delete [] fInputFileNames ; 
266   delete [] fEventNames ; 
267
268   delete fPulse;
269   delete [] fADCValuesLG;
270   delete [] fADCValuesHG;
271
272   if(fcdb){ delete fcdb ; fcdb=0;} 
273
274 }
275
276 //____________________________________________________________________________
277 void AliPHOSDigitizer::Digitize(Int_t event) 
278
279   
280   // Makes the digitization of the collected summable digits.
281   //  It first creates the array of all PHOS modules
282   //  filled with noise (different for EMC, and CPV) and
283   //  then adds contributions from SDigits. 
284   // This design avoids scanning over the list of digits to add 
285   // contribution to new SDigits only.
286
287   Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
288
289   //First stream 
290   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
291   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
292
293   Int_t readEvent = event ; 
294   if (fManager) 
295     readEvent = static_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ; 
296   AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
297                   readEvent, GetTitle(), fEventFolderName.Data())) ; 
298   rl->GetEvent(readEvent) ;
299   phosLoader->CleanSDigits() ; 
300   phosLoader->LoadSDigits("READ") ;
301
302   //Prepare Output
303   TClonesArray * digits = phosLoader->Digits() ;
304   if( !digits ) {
305     phosLoader->MakeDigitsArray() ; 
306     digits = phosLoader->Digits() ;
307   }
308  
309   digits->Clear() ;
310
311   //
312   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
313   //Making digits with noise, first EMC
314   //Check which PHOS modules are present
315   Bool_t isPresent[5] ;
316   TString volpath ;
317   Int_t nmod=0 ;
318   for(Int_t i=0; i<5; i++){
319     volpath = "/ALIC_1/PHOS_";
320     volpath += i+1;
321     if (gGeoManager->CheckPath(volpath.Data())) {
322       isPresent[i]=1 ;
323       nmod++ ;
324     }
325     else{
326       isPresent[i]=0 ;
327     }
328   }
329
330   Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
331   
332   Int_t nCPV ;
333   Int_t absID ;
334   
335   //check if CPV exists
336   Bool_t isCPVpresent=0 ;
337   for(Int_t i=1; i<=5 && !isCPVpresent; i++){
338     volpath = "/ALIC_1/PHOS_";
339     volpath += i;
340     volpath += "/PCPV_1";
341     if (gGeoManager->CheckPath(volpath.Data())) 
342       isCPVpresent=1 ;
343   } 
344   
345   if(isCPVpresent){
346     nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
347   }
348   else{
349      nCPV = nEMC ;
350   }  
351
352   digits->Expand(nCPV) ;
353
354   //take all the inputs to add together and load the SDigits
355   TObjArray * sdigArray = new TObjArray(fInput) ;
356   sdigArray->AddAt(phosLoader->SDigits(), 0) ;
357  
358   for(Int_t i = 1 ; i < fInput ; i++){
359     TString tempo(fEventNames[i]) ; 
360     tempo += i ;
361     AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
362     if(!rl2){ 
363       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
364       if (!rl2) {
365         Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
366         return ;
367       }
368       rl2->LoadHeader();
369     }
370     AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
371  
372     if(fManager){ 
373       readEvent = static_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ; 
374     }
375     TClonesArray * digs ;
376     if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
377       AliInfo(Form("Adding event %d from input stream %d %s %s", 
378                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ; 
379        rl2->GetEvent(readEvent) ;
380        phosLoader2->CleanDigits() ;
381        phosLoader2->LoadDigits("READ") ;
382        digs = phosLoader2->Digits() ;
383        toMakeNoise=0 ; //Do not add noise, it is already in stream
384     }
385     else{
386       AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
387                  readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
388        rl2->GetEvent(readEvent) ;
389        phosLoader2->CleanSDigits() ; 
390        phosLoader2->LoadSDigits("READ") ;
391        digs = phosLoader2->SDigits() ;
392     } 
393     sdigArray->AddAt(digs, i) ;
394   }
395
396   //Find the first crystal with signal
397   Int_t nextSig = 200000 ; 
398   TClonesArray * sdigits ;  
399   for(Int_t i = 0 ; i < fInput ; i++){
400     sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
401     if ( !sdigits->GetEntriesFast() )
402       continue ; 
403     Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
404     if(curNext < nextSig) 
405       nextSig = curNext ;
406   }
407   
408   TArrayI index(fInput) ;
409   index.Reset() ;  //Set all indexes to zero
410   
411   AliPHOSDigit * digit ;
412   AliPHOSDigit * curSDigit ;
413   
414 //  TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
415   
416   //Put Noise contribution
417   Double_t apdNoise = 0. ;
418   if(toMakeNoise)
419      apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ; 
420
421   Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
422   Int_t idigit= 0;
423   for(Int_t imod=0; imod<5; imod++){
424     if(!isPresent[imod])
425       continue ;
426     Int_t firstAbsId=imod*emcpermod+1 ;
427     Int_t lastAbsId =(imod+1)*emcpermod ; 
428     for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
429       Float_t noise = gRandom->Gaus(0.,apdNoise) ; 
430       new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
431       //look if we have to add signal?
432       digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
433       idigit++ ;
434     
435       if(absID==nextSig){
436       //Add SDigits from all inputs 
437 //      ticks->Clear() ;
438 //      Int_t contrib = 0 ;
439
440 //New Timing model is necessary
441 //      Float_t a = digit->GetEnergy() ;
442 //      Float_t b = TMath::Abs( a / fTimeSignalLength) ;
443 //      //Mark the beginning of the signal
444 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);  
445 //      //Mark the end of the signal     
446 //      new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b); 
447
448 // Calculate time as time of the largest digit
449         Float_t time = digit->GetTime() ;
450         Float_t eTime= digit->GetEnergy() ;
451       
452         //loop over inputs
453         for(Int_t i = 0 ; i < fInput ; i++){
454           if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
455             curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;       
456             if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
457               curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
458               curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
459             }
460           }
461           else
462             curSDigit = 0 ;
463           //May be several digits will contribute from the same input
464           while(curSDigit && curSDigit->GetId() == absID){         
465             //Shift primary to separate primaries belonging different inputs
466             Int_t primaryoffset ;
467             if(fManager)
468               primaryoffset = fManager->GetMask(i) ; 
469             else
470               primaryoffset = 10000000*i ;
471             curSDigit->ShiftPrimary(primaryoffset) ;
472           
473 //New Timing model is necessary
474 //        a = curSDigit->GetEnergy() ;
475 //        b = a /fTimeSignalLength ;
476 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);  
477 //        new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
478           if(curSDigit->GetEnergy()>eTime){
479             eTime=curSDigit->GetEnergy() ;
480             time=curSDigit->GetTime() ;
481           }
482             *digit += *curSDigit ;  //add energies
483
484             index[i]++ ;
485             if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
486               curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;     
487             else
488               curSDigit = 0 ;
489           }
490         }
491       
492         //calculate and set time
493 //New Timing model is necessary
494 //      Float_t time = FrontEdgeTime(ticks) ;
495         digit->SetTime(time) ;
496       
497         //Find next signal module
498         nextSig = 200000 ;
499         for(Int_t i = 0 ; i < fInput ; i++){
500           sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
501           Int_t curNext = nextSig ;
502           if(sdigits->GetEntriesFast() > index[i] ){
503             curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
504           }
505           if(curNext < nextSig) nextSig = curNext ;
506         }
507       }
508     }
509   }
510
511   //distretize energy if necessary
512   if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
513     Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
514     for(Int_t i = 0 ; i < nEMC ; i++){                                                                                                       
515       digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
516       digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
517     } 
518   }
519   //Apply decalibration if necessary
520   for(Int_t i = 0 ; i < nEMC ; i++){
521     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
522     Decalibrate(digit) ;
523   }
524  
525   
526 //  ticks->Delete() ;
527 //  delete ticks ;
528   
529   //Now CPV digits (different noise and no timing)
530   Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
531   Int_t nEMCtotal=emcpermod*5 ;
532   Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
533   if(isCPVpresent){  //CPV is present in current geometry
534     for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
535       if(!isPresent[imod])
536         continue ;
537       Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
538       Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
539       for(absID = firstAbsId; absID <= lastAbsId; absID++){
540         Float_t noise = gRandom->Gaus(0., cpvNoise) ; 
541         new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
542         idigit++ ;
543         //look if we have to add signal?
544         if(absID==nextSig){
545           digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
546           //Add SDigits from all inputs
547           for(Int_t i = 0 ; i < fInput ; 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              //May be several digits will contribute from the same input
554              while(curSDigit && curSDigit->GetId() == absID){      
555                //Shift primary to separate primaries belonging different inputs
556                Int_t primaryoffset ;
557                if(fManager)
558                  primaryoffset = fManager->GetMask(i) ; 
559                else
560                  primaryoffset = 10000000*i ;
561                curSDigit->ShiftPrimary(primaryoffset) ;
562
563                //add energies
564                *digit += *curSDigit ;  
565                index[i]++ ;
566                if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
567                  curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;        
568                else
569                  curSDigit = 0 ;
570              }
571           }
572
573           //Find next signal module
574           nextSig = 200000 ;
575           for(Int_t i = 0 ; i < fInput ; i++){
576             sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
577             Int_t curNext = nextSig ;
578             if(sdigits->GetEntriesFast() > index[i] )
579               curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
580             if(curNext < nextSig) nextSig = curNext ;
581           }
582       
583         }
584       }
585     }
586   }
587
588   delete sdigArray ; //We should not delete its contents 
589   
590   Int_t relId[4];
591
592   //set amplitudes in bad channels to zero
593
594   for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
595     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
596     geom->AbsToRelNumbering(digit->GetId(),relId);
597     if(relId[1] == 0) // Emc
598       if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.); 
599   }
600
601   //remove digits below thresholds
602   Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
603   for(Int_t i = 0 ; i < nEMC ; i++){
604     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
605
606     if(digit->GetEnergy() < emcThreshold){
607       digits->RemoveAt(i) ;
608       continue ;
609     }
610
611     geom->AbsToRelNumbering(digit->GetId(),relId);
612
613     digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
614
615     Float_t tres = TimeResolution(digit->GetEnergy()) ; 
616     digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
617
618     fPulse->Reset();
619     fPulse->SetAmplitude(digit->GetEnergy()/
620                          fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
621     fPulse->SetTZero(digit->GetTimeR());
622     fPulse->MakeSamples();
623     fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ; 
624     Int_t nSamples = fPulse->GetRawFormatTimeBins();
625     digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
626     digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
627   }
628
629   Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
630   for(Int_t i = nEMC; i < nCPV ; i++){
631     if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
632       digits->RemoveAt(i) ;
633   } 
634     
635   digits->Compress() ;  
636   Int_t ndigits = digits->GetEntriesFast() ;
637   
638   //Set indexes in list of digits and make true digitization of the energy
639   for (Int_t i = 0 ; i < ndigits ; i++) { 
640     digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; 
641     digit->SetIndexInList(i) ;     
642     if(digit->GetId() > fEmcCrystals){ //digitize CPV only
643       digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
644     }
645   }
646
647 }
648 //____________________________________________________________________________
649 Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
650   //Apply calibration
651   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
652
653   //Determine rel.position of the cell absolute ID
654   Int_t relId[4];
655   geom->AbsToRelNumbering(absId,relId);
656   Int_t module=relId[0];
657   Int_t row   =relId[2];
658   Int_t column=relId[3];
659   if(relId[1]==0){ //This Is EMC
660     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
661     return amp*calibration ;
662   }
663   return 0 ;
664 }
665 //____________________________________________________________________________
666 void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
667 {
668   // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
669
670   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
671
672   //Determine rel.position of the cell absolute ID
673   Int_t relId[4];
674   geom->AbsToRelNumbering(digit->GetId(),relId);
675   Int_t module=relId[0];
676   Int_t row   =relId[2];
677   Int_t column=relId[3];
678   if(relId[1]==0){ //This Is EMC
679     Float_t decalib     = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
680     Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
681     Float_t energy = digit->GetEnergy()/calibration;
682     digit->SetEnergy(energy); //Now digit measures E in ADC counts
683     Float_t time = digit->GetTime() ;
684     time-=fcdb->GetTimeShiftEmc(module,column,row);
685     digit->SetTime(time) ;
686   }
687 }
688 //____________________________________________________________________________
689 Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
690   //Apply time calibration
691   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
692
693   //Determine rel.position of the cell absolute ID
694   Int_t relId[4];
695   geom->AbsToRelNumbering(absId,relId);
696   Int_t module=relId[0];
697   Int_t row   =relId[2];
698   Int_t column=relId[3];
699   time += fcdb->GetTimeShiftEmc(module,column,row);
700   return time ;
701 }
702 //____________________________________________________________________________
703 Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
704 {
705   // Returns digitized value of the CPV charge in a pad absId
706
707   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
708
709   //Determine rel.position of the cell absId
710   Int_t relId[4];
711   geom->AbsToRelNumbering(absId,relId);
712   Int_t module=relId[0];
713   Int_t row   =relId[2];
714   Int_t column=relId[3];
715   
716   Int_t channel = 0;
717   
718   if(absId > fEmcCrystals){ //digitize CPV only
719
720     //reading calibration data for cell absId.
721     Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
722     Float_t adcChanelCpv   = fcdb->GetADCchannelCpv( module,column,row);
723
724     channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;       
725     Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
726     if(channel > nMax ) channel = nMax ;
727   }
728   return channel ;
729 }
730
731 //____________________________________________________________________________
732 void AliPHOSDigitizer::Exec(Option_t *option) 
733
734   // Steering method to process digitization for events
735   // in the range from fFirstEvent to fLastEvent.
736   // This range is optionally set by SetEventRange().
737   // if fLastEvent=-1, then process events until the end.
738   // by default fLastEvent = fFirstEvent (process only one event)
739
740   if (!fInit) { // to prevent overwrite existing file
741     AliError(Form("Give a version name different from %s", 
742                   fEventFolderName.Data() )) ;
743     return ;
744   }   
745
746   if (strstr(option,"print")) {
747     Print();
748     return ; 
749   }
750   
751  if(strstr(option,"tim"))
752      gBenchmark->Start("PHOSDigitizer");
753   
754   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
755   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
756
757   // Post Digitizer to the white board
758   phosLoader->PostDigitizer(this) ;
759   
760   if (fLastEvent == -1) 
761     fLastEvent = rl->GetNumberOfEvents() - 1 ;
762   else if (fManager) 
763     fLastEvent = fFirstEvent ; 
764  
765   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
766   
767   Int_t ievent ;
768
769   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
770     fEventCounter++ ; 
771
772     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
773
774     WriteDigits() ;
775
776     if(strstr(option,"deb"))
777       PrintDigits(option);
778     
779     //increment the total number of Digits per run 
780     fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;  
781  }
782  
783   phosLoader->CleanDigitizer();
784
785   if(strstr(option,"tim")){
786     gBenchmark->Stop("PHOSDigitizer");
787     TString message ; 
788     message = "  took %f seconds for Digitizing %f seconds per event\n" ; 
789     AliInfo(Form( message.Data(), 
790          gBenchmark->GetCpuTime("PHOSDigitizer"), 
791          gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents )); 
792   } 
793 }
794 //____________________________________________________________________________ 
795 Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
796   //calculate TOF resolution using beam-test resutls
797   Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
798   Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
799   return TMath::Sqrt(a*a+b*b/e) ;
800 }
801
802 ////____________________________________________________________________________ 
803 //Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
804 //{
805 //  // Returns the shortest time among all time ticks
806 //
807 //  ticks->Sort() ; //Sort in accordance with times of ticks
808 //  TIter it(ticks) ;
809 //  AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
810 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
811 //
812 //  AliPHOSTick * t ;  
813 //  while((t=(AliPHOSTick*) it.Next())){
814 //    if(t->GetTime() < time)  //This tick starts before crossing
815 //      *ctick+=*t ;
816 //    else
817 //      return time ;
818 //
819 //    time = ctick->CrossingTime(fTimeThreshold) ;    
820 //  }
821 //  return time ;
822 //}
823
824 //____________________________________________________________________________ 
825 Bool_t AliPHOSDigitizer::Init()
826 {
827   // Makes all memory allocations
828   fInit = kTRUE ; 
829
830   AliPHOSGeometry *geom;
831   if (!(geom = AliPHOSGeometry::GetInstance())) 
832         geom = AliPHOSGeometry::GetInstance("IHEP","");
833 //   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
834
835   fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
836   
837   fFirstEvent = 0 ; 
838   fLastEvent = fFirstEvent ; 
839   if (fManager) 
840     fInput = fManager->GetNinputs() ; 
841   else 
842     fInput           = 1 ;
843   
844   fInputFileNames  = new TString[fInput] ; 
845   fEventNames      = new TString[fInput] ; 
846   fInputFileNames[0] = GetTitle() ; 
847   fEventNames[0]     = fEventFolderName.Data() ; 
848   Int_t index ; 
849   for (index = 1 ; index < fInput ; index++) {
850     fInputFileNames[index] = static_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0); 
851     TString tempo = fManager->GetInputFolderName(index) ;
852     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
853   }
854
855   //to prevent cleaning of this object while GetEvent is called
856   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
857   if(!rl){
858     rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ; 
859   }
860   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
861   phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
862
863   return fInit ; 
864 }
865
866 //____________________________________________________________________________ 
867 void AliPHOSDigitizer::InitParameters()
868 {
869   // Set initial parameters Digitizer
870
871   fDigitsInRun  = 0 ; 
872   SetEventRange(0,-1) ;
873   fPulse = new AliPHOSPulseGenerator();
874   fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
875   fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
876     
877 }
878
879 //__________________________________________________________________
880 void AliPHOSDigitizer::Print(const Option_t *)const 
881 {
882   // Print Digitizer's parameters
883   AliInfo(Form("\n------------------- %s -------------", GetName() )) ; 
884   if( strcmp(fEventFolderName.Data(), "") != 0 ){
885     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
886     
887     Int_t nStreams ; 
888     if (fManager) 
889       nStreams =  GetNInputStreams() ;
890     else 
891       nStreams = fInput ; 
892     
893     Int_t index = 0 ;  
894     for (index = 0 ; index < nStreams ; index++) {  
895       TString tempo(fEventNames[index]) ; 
896       tempo += index ;
897       printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ; 
898     }
899     
900  //   printf("\nWith following parameters:\n") ;
901  //   printf("  Electronics noise in EMC (fPinNoise)    = %f GeV\n", fPinNoise ) ; 
902  //   printf("  Threshold  in EMC  (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;  
903  //   printf("  Noise in CPV (fCPVNoise)                = %f aux units\n", fCPVNoise ) ; 
904  //   printf("  Threshold in CPV (fCPVDigitThreshold)   = %f aux units\n",fCPVDigitThreshold ) ;  
905     printf(" ---------------------------------------------------\n") ;   
906   }
907   else
908     AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
909   
910 }
911
912 //__________________________________________________________________
913  void AliPHOSDigitizer::PrintDigits(Option_t * option)
914 {
915   // Print a table of digits
916   
917
918   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
919   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
920   TClonesArray * digits = phosLoader->Digits() ; 
921   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
922   
923   AliInfo(Form("%d", digits->GetEntriesFast())) ; 
924   printf("\nevent %d", gAlice->GetEvNumber()) ;
925   printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;  
926
927  
928   if(strstr(option,"all")||strstr(option,"EMC")){  
929     //loop over digits
930     AliPHOSDigit * digit;
931     printf("\nEMC digits (with primaries):\n")  ;
932     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
933     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
934     Int_t index ;
935     for (index = 0 ; (index < digits->GetEntriesFast()) && 
936            (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
937       digit = (AliPHOSDigit * )  digits->At(index) ;
938       if(digit->GetNprimary() == 0) 
939         continue;
940 //       printf("%6d  %8d    %6.5e %4d      %2d :",
941 //            digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  // YVK
942       printf("%6d  %.4f    %6.5e %4d      %2d :",
943               digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
944       Int_t iprimary;
945       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
946         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
947       }    
948       printf("\n") ;
949     }
950   }
951   
952   if(strstr(option,"all")||strstr(option,"CPV")){
953     
954     //loop over CPV digits
955     AliPHOSDigit * digit;
956     printf("\nCPV digits (with primaries):\n")  ;
957     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
958     Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
959     Int_t index ;
960     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
961       digit = (AliPHOSDigit * )  digits->At(index) ;
962       if(digit->GetId() > maxEmc){
963         printf("%6d  %8d    %4d      %2d :",
964                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;  
965         Int_t iprimary;
966         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
967           printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
968         }    
969         printf("\n") ;
970       }
971     }
972   }
973
974 }
975
976 //__________________________________________________________________
977 Float_t AliPHOSDigitizer::TimeOfNoise(void) const
978 {  // Calculates the time signal generated by noise
979   //PH  Info("TimeOfNoise", "Change me") ; 
980   return gRandom->Rndm() * 1.28E-5;
981 }
982
983 //__________________________________________________________________
984 void AliPHOSDigitizer::Unload() 
985 {  
986   
987   Int_t i ; 
988   for(i = 1 ; i < fInput ; i++){
989     TString tempo(fEventNames[i]) ; 
990     tempo += i ;
991     AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ; 
992     AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
993     phosLoader->UnloadSDigits() ; 
994   }
995   
996   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
997   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
998   phosLoader->UnloadDigits() ; 
999 }
1000
1001 //____________________________________________________________________________
1002 void AliPHOSDigitizer::WriteDigits()
1003 {
1004
1005   // Makes TreeD in the output file. 
1006   // Check if branch already exists: 
1007   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
1008   //      already existing branches. 
1009   //   else creates branch with Digits, named "PHOS", title "...",
1010   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
1011   //      and names of files, from which digits are made.
1012
1013   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1014   AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1015  
1016   const TClonesArray * digits = phosLoader->Digits() ; 
1017   TTree * treeD = phosLoader->TreeD();
1018   if(!treeD){
1019     phosLoader->MakeTree("D");
1020     treeD = phosLoader->TreeD();
1021   }
1022   
1023   // -- create Digits branch
1024   Int_t bufferSize = 32000 ;    
1025   TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1026   digitsBranch->SetTitle(fEventFolderName);
1027   digitsBranch->Fill() ;
1028   
1029   phosLoader->WriteDigits("OVERWRITE");
1030   phosLoader->WriteDigitizer("OVERWRITE");
1031
1032   Unload() ; 
1033
1034 }