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