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