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