Added a protection in the dtor. When the tasks is created by default ctor (to access...
[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 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute) 
20 //////////////////////////////////////////////////////////////////////////////
21 // This TTask performs digitization of Summable digits (in the PHOS case it is just
22 // the sum of contributions from all primary particles into a given cell). 
23 // In addition it performs mixing of summable digits from different events.
24 // The name of the TTask is also the title of the branch that will contain 
25 // the created SDigits
26 // The title of the TTAsk is the name of the file that contains the hits from
27 // which the SDigits are created
28 //
29 // For each event two branches are created in TreeD:
30 //   "PHOS" - list of digits
31 //   "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
32 //
33 // Note, that one can set a title for new digits branch, and repeat digitization with
34 // another set of parameters.
35 //
36 // Use case:
37 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
38 // root[1] d->ExecuteTask()             
39 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
40 //                       //Digitizes SDigitis in all events found in file galice.root 
41 //
42 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;  
43 //                       // Will read sdigits from galice1.root
44 // root[3] d1->MixWith("galice2.root")       
45 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
46 //                       // Reads another set of sdigits from galice2.root
47 // root[3] d1->MixWith("galice3.root")       
48 //                       // Reads another set of sdigits from galice3.root
49 // root[4] d->ExecuteTask("deb timing")    
50 //                       // Reads SDigits from files galice1.root, galice2.root ....
51 //                       // mixes them and stores produced Digits in file galice1.root          
52 //                       // deb - prints number of produced digits
53 //                       // deb all - prints list of produced digits
54 //                       // timing  - prints time used for digitization
55 //
56
57 // --- ROOT system ---
58 #include "TFile.h"
59 #include "TTree.h"
60 #include "TSystem.h"
61 #include "TROOT.h"
62 #include "TFolder.h"
63 #include "TObjString.h"
64 #include "TGeometry.h"
65 #include "TBenchmark.h"
66
67 // --- Standard library ---
68 #include <iomanip.h>
69
70 // --- AliRoot header files ---
71
72 #include "AliRun.h"
73 #include "AliHeader.h"
74 #include "AliRunDigitizer.h"
75 #include "AliPHOSDigit.h"
76 #include "AliPHOS.h"
77 #include "AliPHOSGetter.h"
78 #include "AliPHOSDigitizer.h"
79 #include "AliPHOSSDigitizer.h"
80 #include "AliPHOSGeometry.h"
81 #include "AliPHOSTick.h"
82
83 ClassImp(AliPHOSDigitizer)
84
85
86 //____________________________________________________________________________ 
87   AliPHOSDigitizer::AliPHOSDigitizer() 
88 {
89   // ctor
90
91   InitParameters() ; 
92 }
93
94 //____________________________________________________________________________ 
95 AliPHOSDigitizer::AliPHOSDigitizer(const char *headerFile,const char * name)
96 {
97   // ctor
98   SetTitle(headerFile) ;
99   SetName(name) ;
100   fManager = 0 ;                     // We work in the standalong mode
101   fSplitFile= 0 ; 
102   InitParameters() ; 
103   Init() ;
104   
105 }
106
107 //____________________________________________________________________________ 
108 AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * ard):AliDigitizer(ard)
109 {
110   // ctor
111   SetTitle("aliroot") ;
112   SetName("Default") ;
113   InitParameters() ; 
114   
115 }
116
117 //____________________________________________________________________________ 
118   AliPHOSDigitizer::~AliPHOSDigitizer()
119 {
120   // dtor
121   // gime=0 if Digitizer created by default ctor (to get just the parameters)
122   
123  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
124  
125  if (gime) {
126    // remove the task from the folder list
127    gime->RemoveTask("S",GetName()) ;
128    gime->RemoveTask("D",GetName()) ;
129    
130    // remove the Digits from the folder list
131    gime->RemoveObjects("D", GetName()) ;
132    
133    // remove the SDigits from the folder list
134    gime->RemoveSDigits() ;
135    
136    // Delete gAlice
137    gime->CloseFile() ; 
138    
139    fSplitFile = 0 ; 
140  }
141 }
142
143 //____________________________________________________________________________
144 void AliPHOSDigitizer::Digitize(const Int_t event) 
145
146   
147   // Makes the digitization of the collected summable digits.
148   //  It first creates the array of all PHOS modules
149   //  filled with noise (different for EMC, CPV and PPSD) and
150   //  then adds contributions from SDigits. 
151   // This design avoids scanning over the list of digits to add 
152   // contribution to new SDigits only.
153
154   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
155   TClonesArray * digits = gime->Digits(GetName()) ; 
156
157   digits->Clear() ;
158
159   const AliPHOSGeometry *geom = gime->PHOSGeometry() ; 
160   //Making digits with noise, first EMC
161   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
162   
163   Int_t nCPV ;
164   Int_t absID ;
165   TString name      =  geom->GetName() ;
166   
167   nCPV = nEMC + geom->GetNumberOfCPVPadsZ()*geom->GetNumberOfCPVPadsPhi()*
168     geom->GetNModules() ;
169
170   digits->Expand(nCPV) ;
171
172   // get first the sdigitizer from the tasks list (must have same name as the digitizer)
173   const AliPHOSSDigitizer * sDigitizer = new AliPHOSSDigitizer() ; //gime->SDigitizer(GetName()); 
174   if ( !sDigitizer) {
175     cerr << "ERROR: AliPHOSDigitizer::Digitize -> SDigitizer with name " << GetName() << " not found " << endl ; 
176     abort() ; 
177   }
178     
179   // loop through the sdigits posted to the White Board and add them to the noise
180   TCollection * folderslist = gime->SDigitsFolder()->GetListOfFolders() ; 
181   TIter next(folderslist) ; 
182   TFolder * folder = 0 ; 
183   TClonesArray * sdigits = 0 ;
184   Int_t input = 0 ;
185   TObjArray * sdigArray = new TObjArray(2) ;
186   while ( (folder = (TFolder*)next()) ) { 
187     if ( (sdigits = (TClonesArray*)folder->FindObject(GetName()) ) ) {
188       TString fileName(folder->GetName()) ;
189       fileName.ReplaceAll("_","/") ;
190       cout << "INFO: AliPHOSDigitizer::Digitize -> Adding SDigits " 
191            << GetName() << " from " << fileName << endl ; 
192       sdigArray->AddAt(sdigits, input) ;
193       input++ ;
194     }
195   }
196
197   //Find the first crystall with signal
198   Int_t nextSig = 200000 ; 
199   Int_t i;
200   for(i=0; i<input; i++){
201     sdigits = (TClonesArray *)sdigArray->At(i) ;
202     if ( !sdigits->GetEntriesFast() )
203       continue ; 
204     Int_t curNext = ((AliPHOSDigit *)sdigits->At(0))->GetId() ;
205      if(curNext < nextSig) 
206        nextSig = curNext ;
207   }
208
209   TArrayI index(input) ;
210   index.Reset() ;  //Set all indexes to zero
211
212   AliPHOSDigit * digit ;
213   AliPHOSDigit * curSDigit ;
214
215   TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
216
217   //Put Noise contribution
218   for(absID = 1; absID <= nEMC; absID++){
219     Float_t noise = gRandom->Gaus(0., fPinNoise) ; 
220     new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
221     //look if we have to add signal?
222     digit = (AliPHOSDigit *) digits->At(absID-1) ;
223  
224     if(absID==nextSig){
225       //Add SDigits from all inputs 
226       ticks->Clear() ;
227       Int_t contrib = 0 ;
228       Float_t a = digit->GetAmp() ;
229       Float_t b = TMath::Abs( a /fTimeSignalLength) ;
230       //Mark the beginnign of the signal
231       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);  
232       //Mark the end of the ignal     
233       new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b); 
234
235       //loop over inputs
236       for(i=0; i<input; i++){
237         if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
238           curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;         
239         else
240           curSDigit = 0 ;
241         //May be several digits will contribute from the same input
242         while(curSDigit && curSDigit->GetId() == absID){           
243           //Shift primary to separate primaries belonging different inputs
244           Int_t primaryoffset ;
245           if(fManager)
246             primaryoffset = fManager->GetMask(i) ; 
247           else
248             primaryoffset = 10000000*i ;
249           curSDigit->ShiftPrimary(primaryoffset) ;
250           
251           a = curSDigit->GetAmp() ;
252           b = a /fTimeSignalLength ;
253           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);  
254           new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
255
256           *digit = *digit + *curSDigit ;  //add energies
257
258           index[i]++ ;
259           if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
260             curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;       
261           else
262             curSDigit = 0 ;
263         }
264       }
265
266       //calculate and set time
267       Float_t time = FrontEdgeTime(ticks) ;
268       digit->SetTime(time) ;
269
270       //Find next signal module
271       nextSig = 200000 ;
272       for(i=0; i<input; i++){
273         sdigits = ((TClonesArray *)sdigArray->At(i)) ;
274         Int_t curNext = nextSig ;
275         if(sdigits->GetEntriesFast() > index[i] ){
276           curNext = ((AliPHOSDigit *) sdigits->At(index[i]))->GetId() ;
277           
278         }
279         if(curNext < nextSig) nextSig = curNext ;
280       }
281     }
282   }
283   
284   ticks->Delete() ;
285   delete ticks ;
286
287
288   //Now CPV digits (different noise and no timing)
289   for(absID = nEMC+1; absID <= nCPV; absID++){
290     Float_t noise = gRandom->Gaus(0., fCPVNoise) ; 
291     new((*digits)[absID-1]) AliPHOSDigit( -1,absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
292     //look if we have to add signal?
293     if(absID==nextSig){
294       digit = (AliPHOSDigit *) digits->At(absID-1) ;
295       //Add SDigits from all inputs
296       for(i=0; i<input; i++){
297         if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
298           curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;         
299         else
300           curSDigit = 0 ;
301
302         //May be several digits will contribute from the same input
303         while(curSDigit && curSDigit->GetId() == absID){           
304           //Shift primary to separate primaries belonging different inputs
305           Int_t primaryoffset ;
306           if(fManager)
307             primaryoffset = fManager->GetMask(i) ; 
308           else
309             primaryoffset = 10000000*i ;
310           curSDigit->ShiftPrimary(primaryoffset) ;
311
312           //add energies
313           *digit = *digit + *curSDigit ;  
314           index[i]++ ;
315           if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
316             curSDigit = (AliPHOSDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;       
317           else
318             curSDigit = 0 ;
319         }
320       }
321
322       //Find next signal module
323       nextSig = 200000 ;
324       for(i=0; i<input; i++){
325         sdigits = (TClonesArray *)sdigArray->At(i) ;
326         Int_t curNext = nextSig ;
327         if(sdigits->GetEntriesFast() > index[i] )
328           curNext = ((AliPHOSDigit *) sdigits->At(index[i]))->GetId() ;
329         if(curNext < nextSig) nextSig = curNext ;
330       }
331       
332     }
333   }
334   delete sdigArray ; //We should not delete its contents
335   
336   
337   
338   //remove digits below thresholds
339   for(i = 0; i < nEMC ; i++){
340     digit = (AliPHOSDigit*) digits->At(i) ;
341     if(sDigitizer->Calibrate( digit->GetAmp() ) < fEMCDigitThreshold)
342       digits->RemoveAt(i) ;
343     else
344       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
345   }
346
347
348   for(i = nEMC; i < nCPV ; i++)
349     if(sDigitizer->Calibrate(((AliPHOSDigit*)digits->At(i))->GetAmp()) < fCPVDigitThreshold)
350       digits->RemoveAt(i) ;
351     
352   digits->Compress() ;  
353   
354   Int_t ndigits = digits->GetEntriesFast() ;
355   digits->Expand(ndigits) ;
356
357   //Set indexes in list of digits and make true digitization of the energy
358   for (i = 0 ; i < ndigits ; i++) { 
359     AliPHOSDigit * digit = (AliPHOSDigit *) digits->At(i) ; 
360     digit->SetIndexInList(i) ;     
361     Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ;
362     digit->SetAmp(DigitizeEnergy(energy,digit->GetId()) ) ;
363   }
364
365 }
366 //____________________________________________________________________________
367 Int_t AliPHOSDigitizer::DigitizeEnergy(Float_t energy, Int_t absId)
368 {
369   Int_t chanel ;
370   if(absId <= fEmcCrystals){ //digitize as EMC 
371     chanel = (Int_t) TMath::Ceil((energy - fADCpedestalEmc)/fADCchanelEmc) ;       
372     if(chanel > fNADCemc ) chanel =  fNADCemc ;
373   }
374   else{ //Digitize as CPV
375     chanel = (Int_t) TMath::Ceil((energy - fADCpedestalCpv)/fADCchanelCpv) ;       
376     if(chanel > fNADCcpv ) chanel =  fNADCcpv ;
377   }
378   return chanel ;
379 }
380 //____________________________________________________________________________
381 void AliPHOSDigitizer::Exec(Option_t *option) 
382
383   // Managing method
384
385   if(strcmp(GetName(), "") == 0 )   
386     Init() ;
387   if (strstr(option,"print")) {
388     Print("");
389     return ; 
390   }
391   
392   if(strstr(option,"tim"))
393     gBenchmark->Start("PHOSDigitizer");
394
395   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
396   
397   Int_t nevents ;
398   
399   TTree * treeD ;
400   
401   if(fManager){
402     treeD = fManager->GetTreeD() ;
403     nevents = 1 ;    // Will process only one event
404   }
405   else {
406     gAlice->GetEvent(0) ;
407     nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
408     treeD=gAlice->TreeD() ;
409   }
410
411
412   //Check, if this branch already exits
413   if (treeD) { 
414     TObjArray * lob = (TObjArray*)treeD->GetListOfBranches() ;
415     TIter next(lob) ; 
416     TBranch * branch = 0 ;  
417     Bool_t phosfound = kFALSE, digitizerfound = kFALSE ; 
418     
419     while ( (branch = (TBranch*)next()) && (!phosfound || !digitizerfound) ) {
420       if ( (strcmp(branch->GetName(), "PHOS")==0) && 
421            (strcmp(branch->GetTitle(), GetName())==0) ) 
422         phosfound = kTRUE ;
423       
424       else if ( (strcmp(branch->GetName(), "AliPHOSDigitizer")==0) && 
425                 (strcmp(branch->GetTitle(), GetName())==0) ) 
426         digitizerfound = kTRUE ; 
427     }
428     
429     if ( phosfound ) {
430       cerr << "WARNING: AliPHOSDigitizer -> Digits branch with name " << GetName() 
431            << " already exits" << endl ;
432       return ; 
433     }   
434     if ( digitizerfound ) {
435       cerr << "WARNING: AliPHOSDigitizer -> Digitizer branch with name " << GetName() 
436            << " already exits" << endl ;
437       return ; 
438     }   
439   }
440
441   Int_t ievent ;
442   
443   for(ievent = 0; ievent < nevents; ievent++){
444   
445     if(fManager){
446
447       Int_t input ;
448       for(input = 0 ; input < fManager->GetNinputs(); input ++){
449         TTree * treeS = fManager->GetInputTreeS(input) ;
450         if(!treeS){
451           cerr << "AliPHOSDigitizer::Exec -> No Input " << endl ;
452           return ;
453         }
454         gime->ReadTreeS(treeS,input) ;
455       }
456
457     }
458     else
459       gime->Event(ievent,"S") ;
460     
461     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
462     
463     WriteDigits(ievent) ;
464     
465     if(strstr(option,"deb"))
466       PrintDigits(option);
467     
468     //increment the total number of Digits per run 
469     fDigitsInRun += gime->Digits()->GetEntriesFast() ;  
470   }
471    
472   if(strstr(option,"tim")){
473     gBenchmark->Stop("PHOSDigitizer");
474     cout << "AliPHOSDigitizer:" << endl ;
475     cout << "  took " << gBenchmark->GetCpuTime("PHOSDigitizer") << " seconds for Digitizing " 
476          <<  gBenchmark->GetCpuTime("PHOSDigitizer")/nevents << " seconds per event " << endl ;
477     cout << endl ;
478   }
479   
480 }
481
482 //____________________________________________________________________________ 
483 Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) 
484 { // 
485   ticks->Sort() ; //Sort in accordance with times of ticks
486   TIter it(ticks) ;
487   AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
488   Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
489
490   AliPHOSTick * t ;  
491   while((t=(AliPHOSTick*) it.Next())){
492     if(t->GetTime() < time)  //This tick starts before crossing
493       *ctick+=*t ;
494     else
495       return time ;
496
497     time = ctick->CrossingTime(fTimeThreshold) ;    
498   }
499   return time ;
500 }
501
502 //____________________________________________________________________________ 
503 Bool_t AliPHOSDigitizer::Init()
504 {
505
506   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), GetName()) ; 
507   if ( gime == 0 ) {
508     cerr << "ERROR: AliPHOSDigitizer::Init -> Could not obtain the Getter object !" << endl ; 
509     return kFALSE;
510   } 
511   
512   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
513
514   fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
515   
516   // Post Digits to the white board
517   gime->PostDigits(GetName() ) ;   
518   
519   // Post Digitizer to the white board
520   gime->PostDigitizer(this) ;
521   
522   //Mark that we will use current header file
523   if(!fManager){
524     gime->PostSDigits(GetName(),GetTitle()) ;
525     gime->PostSDigitizer(GetName(),GetTitle()) ;
526   }
527   return kTRUE ;
528 }
529
530 //____________________________________________________________________________ 
531 void AliPHOSDigitizer::InitParameters()
532 {
533   fPinNoise           = 0.004 ;
534   fEMCDigitThreshold  = 0.012 ;
535   fCPVNoise           = 0.01;
536   fCPVDigitThreshold  = 0.09 ;
537   fTimeResolution     = 0.5e-9 ;
538   fTimeSignalLength   = 1.0e-9 ;
539   fDigitsInRun  = 0 ; 
540   fADCchanelEmc = 0.0015;        // width of one ADC channel in GeV
541   fADCpedestalEmc = 0.005 ;      //
542   fNADCemc = (Int_t) TMath::Power(2,16) ;  // number of channels in EMC ADC
543
544   fADCchanelCpv = 0.0012 ;          // width of one ADC channel in CPV 'popugais'
545   fADCpedestalCpv = 0.012 ;         // 
546   fNADCcpv = (Int_t) TMath::Power(2,12);      // number of channels in CPV ADC
547
548   fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
549     
550 }
551
552 //__________________________________________________________________
553 void AliPHOSDigitizer::MixWith(const char* headerFile)
554 {
555   // Allows to produce digits by superimposing background and signal event.
556   // It is assumed, that headers file with SIGNAL events is opened in 
557   // the constructor. 
558   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed 
559   // Thus we avoid writing (changing) huge and expensive 
560   // backgound files: all output will be writen into SIGNAL, i.e. 
561   // opened in constructor file. 
562   //
563   // One can open as many files to mix with as one needs.
564   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
565   // can be mixed.
566
567   if( strcmp(GetName(), "") == 0 )
568     Init() ;
569
570   if(fManager){
571     cout << "Can not use this method under AliRunDigitizer " << endl ;
572     return ;
573   }
574   
575   // check if the specified SDigits do not already exist on the White Board:
576   // //Folders/RunMC/Event/Data/PHOS/SDigits/headerFile/sdigitsname
577
578   TString path = "Folders/RunMC/Event/Data/PHOS/SDigits" ; 
579   path += headerFile ; 
580   path += "/" ; 
581   path += GetName() ;
582   if ( gROOT->FindObjectAny(path.Data()) ) {
583     cerr << "WARNING: AliPHOSDigitizer::MixWith -> Entry already exists, do not add" << endl ;
584     return;
585   }
586
587   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
588   gime->PostSDigits(GetName(),headerFile) ;
589   
590   // check if the requested file is already open or exist and if SDigits Branch exist
591   TFile * file = (TFile*)gROOT->FindObject(headerFile); 
592   if ( !file ) { 
593     file = new TFile(headerFile, "READ") ; 
594     if (!file) { 
595       cerr << "ERROR: AliPHOSDigitizer::MixWith -> File " << headerFile << " does not exist!" << endl ; 
596       return ; 
597     }
598   }
599   
600 }
601
602 //__________________________________________________________________
603 void AliPHOSDigitizer::SetSplitFile(const TString splitFileName)
604 {
605   // Diverts the Digits in a file separate from the hits file
606   
607   // I guess it is not going to work if we do merging
608   if (fManager) {
609     cerr << "ERROR: AliPHOSDigitizer::SetSplitFile -> Not yet available in case of merging activated " << endl ;  
610     return ; 
611   }
612
613   TDirectory * cwd = gDirectory ;
614   if ( !(gAlice->GetTreeDFileName() == splitFileName) ) {
615     if (gAlice->GetTreeDFile() )  
616       gAlice->GetTreeDFile()->Close() ; 
617   }
618
619   fSplitFile = gAlice->InitTreeFile("D",splitFileName.Data());
620   fSplitFile->cd() ; 
621   gAlice->Write(0, TObject::kOverwrite);
622
623   TTree *treeE  = gAlice->TreeE();
624   if (!treeE) {
625     cerr << "ERROR: AliPHOSDigitizer::SetSplitFile -> No TreeE found "<<endl;
626     abort() ;
627   }      
628   
629   // copy TreeE
630   AliHeader *header = new AliHeader();
631   treeE->SetBranchAddress("Header", &header);
632   treeE->SetBranchStatus("*",1);
633   TTree *treeENew =  treeE->CloneTree();
634   treeENew->Write(0, TObject::kOverwrite);
635   
636   // copy AliceGeom
637   TGeometry *AliceGeom = static_cast<TGeometry*>(cwd->Get("AliceGeom"));
638   if (!AliceGeom) {
639     cerr << "ERROR: AliPHOSDigitizer::SetSplitFile -> AliceGeom was not found in the input file "<<endl;
640     abort() ;
641     }
642   AliceGeom->Write(0, TObject::kOverwrite);
643   
644   gAlice->MakeTree("D",fSplitFile);
645   cwd->cd() ; 
646   cout << "INFO: AliPHOSDigitizer::SetSPlitMode -> Digits will be stored in " << splitFileName.Data() << endl ; 
647 }
648
649 //__________________________________________________________________
650 void AliPHOSDigitizer::Print(Option_t* option)const {
651   // Print Digitizer's parameters
652   if( strcmp(GetName(), "") != 0 ){
653     
654     cout << "------------------- "<< GetName() << " -------------" << endl ;
655     cout << "Digitizing sDigits from file(s): " <<endl ;
656     
657      TCollection * folderslist = ((TFolder*)gROOT->FindObjectAny("Folders/RunMC/Event/Data/PHOS/SDigits"))->GetListOfFolders() ; 
658     TIter next(folderslist) ; 
659     TFolder * folder = 0 ; 
660     
661     while ( (folder = (TFolder*)next()) ) {
662       if ( folder->FindObject(GetName())  ) 
663         cout << "Adding SDigits " << GetName() << " from " << folder->GetName() << endl ; 
664     }
665     cout << endl ;
666     cout << "Writing digits to " << GetTitle() << endl ;
667     
668     cout << endl ;
669     cout << "With following parameters: " << endl ;
670     cout << "     Electronics noise in EMC (fPinNoise) = " << fPinNoise << endl ;
671     cout << "  Threshold  in EMC  (fEMCDigitThreshold) = " << fEMCDigitThreshold  << endl ; ;
672     cout << "                 Noise in CPV (fCPVNoise) = " << fCPVNoise << endl ; 
673     cout << "    Threshold in CPV (fCPVDigitThreshold) = " << fCPVDigitThreshold << endl ; 
674     cout << "---------------------------------------------------" << endl ;
675   }
676   else
677     cout << "AliPHOSDigitizer not initialized " << endl ;
678   
679 }
680
681 //__________________________________________________________________
682 void AliPHOSDigitizer::PrintDigits(Option_t * option){
683   // Print a table of digits
684
685   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
686   TClonesArray * digits = gime->Digits() ; 
687
688   cout << "AliPHOSDigitiser: event " << gAlice->GetEvNumber() << endl ;
689   cout << "       Number of entries in Digits list " << digits->GetEntriesFast() << endl ;
690   cout << endl ;
691   if(strstr(option,"all")||strstr(option,"EMC")){
692     
693     //loop over digits
694     AliPHOSDigit * digit;
695     cout << "EMC digits (with primaries): " << endl ;
696     cout << "Digit Id     Amplitude   Index       Nprim  Primaries list " <<  endl;      
697     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
698     Int_t index ;
699     for (index = 0 ; (index < digits->GetEntriesFast()) && 
700          (((AliPHOSDigit * )  digits->At(index))->GetId() <= maxEmc) ; index++) {
701       digit = (AliPHOSDigit * )  digits->At(index) ;
702       if(digit->GetNprimary() == 0) continue;
703       cout << setw(6)  <<  digit->GetId() << "   "  <<  setw(10)  <<  digit->GetAmp() <<   "    "  
704            << setw(6)  <<  digit->GetIndexInList() << "    "   
705            << setw(5)  <<  digit->GetNprimary() <<"    ";
706       
707       Int_t iprimary;
708       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
709         cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "    ";
710       cout << endl;      
711     }    
712     cout << endl;
713   }
714
715   if(strstr(option,"all")||strstr(option,"CPV")){
716     
717     //loop over CPV digits
718     AliPHOSDigit * digit;
719     cout << "CPV digits: " << endl ;
720     cout << "Digit Id     Amplitude   Index       Nprim  Primaries list " <<  endl;      
721     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
722     Int_t index ;
723     for (index = 0 ; index < digits->GetEntriesFast(); index++) {
724       digit = (AliPHOSDigit * )  digits->At(index) ;
725       if(digit->GetId() > maxEmc){
726         cout << setw(6)  <<  digit->GetId() << "   "  <<        setw(10)  <<  digit->GetAmp() <<   "    "  
727              << setw(6)  <<  digit->GetIndexInList() << "    "   
728              << setw(5)  <<  digit->GetNprimary() <<"    ";
729         
730         Int_t iprimary;
731         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
732           cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "    ";
733         cout << endl;    
734       }    
735     }
736   }
737
738 }
739
740 //__________________________________________________________________
741 void AliPHOSDigitizer::SetSDigitsBranch(const char* title)
742 {
743   // we set title (comment) of the SDigits branch in the first! header file
744   if( strcmp(GetName(), "") == 0 )
745     Init() ;
746
747   AliPHOSGetter::GetInstance()->SDigits()->SetName(title) ; 
748  
749 }
750 //__________________________________________________________________
751 Float_t AliPHOSDigitizer::TimeOfNoise(void)
752 {  // Calculates the time signal generated by noise
753   //to be rewritten, now returns just big number
754   return 1. ;
755
756 }
757 //____________________________________________________________________________
758 void AliPHOSDigitizer::Reset() 
759
760   // sets current event number to the first simulated event
761
762   if( strcmp(GetName(), "") == 0 )
763     Init() ;
764
765  //  Int_t inputs ;
766 //   for(inputs = 0; inputs < fNinputs ;inputs++)
767 //       fIevent->AddAt(-1, inputs ) ;
768   
769 }
770
771 //____________________________________________________________________________
772 void AliPHOSDigitizer::WriteDigits(Int_t event)
773 {
774
775   // Makes TreeD in the output file. 
776   // Check if branch already exists: 
777   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
778   //      already existing branches. 
779   //   else creates branch with Digits, named "PHOS", title "...",
780   //      and branch "AliPHOSDigitizer", with the same title to keep all the parameters
781   //      and names of files, from which digits are made.
782
783   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
784   const TClonesArray * digits = gime->Digits(GetName()) ; 
785   TTree * treeD ;
786
787  
788  if(fManager) 
789    treeD = fManager->GetTreeD() ;
790  else {
791    if (!gAlice->TreeD() ) 
792      gAlice->MakeTree("D", fSplitFile);
793    treeD = gAlice->TreeD();
794  }
795
796
797   // -- create Digits branch
798   Int_t bufferSize = 32000 ;    
799   TBranch * digitsBranch = treeD->Branch("PHOS",&digits,bufferSize);
800   digitsBranch->SetTitle(GetName());
801     
802   // -- Create Digitizer branch
803   Int_t splitlevel = 0 ;
804   AliPHOSDigitizer * d = gime->Digitizer(GetName()) ;
805   TBranch * digitizerBranch = treeD->Branch("AliPHOSDigitizer", "AliPHOSDigitizer", &d,bufferSize,splitlevel); 
806   digitizerBranch->SetTitle(GetName());
807
808   digitsBranch->Fill() ;
809   digitizerBranch->Fill() ; 
810   treeD->AutoSave() ; 
811  
812 }