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