new classes for digitization derived from TTask
[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 // This is a TTask that constructs SDigits out of Hits
20 // A Summable Digits is the sum of all hits in a cell
21 // A threshold is applied 
22 //
23 //*-- Author :  Dmitri Peressounko (SUBATECH & Kurchatov Institute) 
24 //////////////////////////////////////////////////////////////////////////////
25 // Class performs digitization of Summable digits (in the PHOS case this is just
26 // sum of contributions of all primary particles into give cell). 
27 // In addition it performs mixing of summable digits from different events.
28 // Examples of use:
29 // root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
30 // root[1] d->ExecuteTask()             
31 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
32 //                       //Digitizes SDigitis in all events found in file galice.root 
33 //                       //Depending on variable "CONFIG_SPLIT_FILE" reads branches stored in galice.root
34 //                       //or in PHOS.SDigits.root
35 // root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;  // Will read sdigits from galice1.root
36 // root[3] d1->MixWith("galice2.root",1)       // Reads another portion of sdigits from galice2.root
37 //                                             // says, that this will be output file
38 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
39 // root[3] d1->MixWith("galice3.root",1)       // Reads another portion of sdigits from galice3.root
40 //                                             // overwrides previous definition of output file
41 // root[4] d->ExecuteTask()    // Reads SDigits from files galice1.root, galice2.root ....
42 //                             // mixes them and stores produced Digits in file galice3.root          
43 //
44 //
45 // 
46
47 // --- ROOT system ---
48 #include "TTask.h"
49 #include "TTree.h"
50 #include "TSystem.h"
51 // --- Standard library ---
52
53 // --- AliRoot header files ---
54
55 #include "AliPHOSDigit.h"
56 #include "AliPHOSHit.h"
57 #include "AliPHOSv1.h"
58 #include "AliPHOSDigitizer.h"
59 #include "AliPHOSSDigitizer.h"
60 #include "TROOT.h"
61 #include "TFolder.h"
62
63 ClassImp(AliPHOSDigitizer)
64
65
66 //____________________________________________________________________________ 
67   AliPHOSDigitizer::AliPHOSDigitizer():TTask("AliPHOSDigitizer","") 
68 {
69   // ctor
70
71   fSDigitizer = 0 ;
72   fNinputs = 1 ;
73   fPinNoise = 0.01 ;
74   fEMCDigitThreshold = 0.01 ;
75   fCPVNoise = 0.01;
76   fCPVDigitThreshold = 0.09 ;
77   fPPSDNoise = 0.0000001;
78   fPPSDDigitThreshold = 0.0000002 ;
79   fInitialized = kFALSE ;
80   // add Task to //root/Tasks folder
81   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
82   roottasks->Add(this) ; 
83
84 }
85 //____________________________________________________________________________ 
86 void AliPHOSDigitizer::Init(Int_t isOutFile){
87 // Mades all memory allocations and defiles, 
88 // whether first (default) file will be output file (isOutFile !=0) 
89
90   if(!fInitialized){
91
92     fHeaderFiles  = new TClonesArray("TObjString",1) ;
93     new((*fHeaderFiles)[0]) TObjString("galice.root") ;
94     TFile * file ;
95
96     if(isOutFile)
97       file = new TFile(((TObjString *) fHeaderFiles->At(0))->GetString(),"update") ;
98     else
99       file = new TFile(((TObjString *) fHeaderFiles->At(0))->GetString()) ;
100
101     file->cd() ;
102   
103     fSDigitsFiles = new TClonesArray("TObjString",1);
104     if(gSystem->Getenv("CONFIG_SPLIT_FILE")) 
105       new((*fSDigitsFiles)[0]) TObjString("./PHOS.SDigits.root") ;   
106     else
107       new((*fSDigitsFiles)[0]) TObjString("") ;   
108     
109     fSDigits      = new TClonesArray("TClonesArray",1) ;
110     new((*fSDigits)[0]) TClonesArray("AliPHOSDigit",1000) ;
111     
112     fDigits = new TClonesArray("AliPHOSDigit",200000) ;
113     
114     fIevent    = new TArrayI(1) ;
115     fIevent->AddAt(-1,0 ) ; 
116     fIeventMax = new TArrayI(1) ;
117
118     //Store digits in this file
119     if(isOutFile){
120       gAlice = (AliRun *) file->Get("gAlice") ;
121       fIeventMax->AddAt((Int_t) gAlice->TreeE()->GetEntries(), 0 );
122       fOutFileNumber = 0 ;
123     }
124     else{
125       TTree * te = (TTree *) file->Get("TE") ;
126       fIeventMax->AddAt((Int_t) te->GetEntries(), 0 );
127       fOutFileNumber = -1 ;
128     }
129
130     fInitialized = kTRUE ;
131   }
132
133 }
134
135
136 //____________________________________________________________________________ 
137 AliPHOSDigitizer::AliPHOSDigitizer(char *HeaderFile,char *DigitsFile = 0):TTask("AliPHOSDigitizer","")
138 {
139   // ctor
140   fHeaderFiles  = new TClonesArray("TFile",1) ;          
141   new((*fHeaderFiles)[0]) TObjString(HeaderFile) ;
142   TFile * file = new TFile(((TObjString *) fHeaderFiles->At(0))->GetString(),"update") ;      // Header file, where result will be stored
143
144   file->cd() ;
145   
146   fSDigitsFiles = new TClonesArray("TObjString",1);         // File name of the SDigits branch
147   if(DigitsFile)
148     new((*fSDigitsFiles)[0]) TObjString(DigitsFile) ;   
149   else
150     if(gSystem->Getenv("CONFIG_SPLIT_FILE")) 
151       new((*fSDigitsFiles)[0]) TObjString("./PHOS.SDigits.root") ;   
152     else
153       new((*fSDigitsFiles)[0]) TObjString("") ;   
154     
155   fSDigits      = new TClonesArray("TClonesArray",1) ;      // here list of SDigits wil be stored
156   new((*fSDigits)[0]) TClonesArray("AliPHOSDigit",1000) ;
157     
158   fDigits = new TClonesArray("AliPHOSDigit",200000) ;
159   fDigitsFile="PHOS.Digits" ; 
160     
161   fIevent    = new TArrayI(1) ;
162   fIevent->AddAt(-1,0 ) ; 
163   fIeventMax = new TArrayI(1) ;
164   //Should be check whether gAlice in memory is the same as in file
165   //However, there is no such method (?) ==> we are forced to read it
166   // if(gAlice->TreeE()==0)
167
168   gAlice = (AliRun *) file->Get("gAlice") ;  //If not read yet
169
170   // Get number of events to process
171   fIeventMax->AddAt((Int_t) gAlice->TreeE()->GetEntries(), 0 );
172   fOutFileNumber = 0 ;
173
174   fNinputs = 1 ;
175
176   fPinNoise = 0.01 ;
177   fEMCDigitThreshold = 0.01 ;
178   fCPVNoise = 0.01;
179   fCPVDigitThreshold = 0.09 ;
180   fPPSDNoise = 0.0000001;
181   fPPSDDigitThreshold = 0.0000002 ;  
182
183   // add Task to //root/Tasks folder
184   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
185   roottasks->Add(this) ; 
186   fInitialized = kTRUE ;
187     
188 }
189
190 //____________________________________________________________________________ 
191   AliPHOSDigitizer::~AliPHOSDigitizer()
192 {
193   // dtor
194   delete fHeaderFiles ;
195   delete fSDigitsFiles ;
196   delete fSDigits ;
197   delete fDigits ;
198 }
199 //____________________________________________________________________________
200 Bool_t AliPHOSDigitizer::Combinator() { 
201
202   //Makes all desirable combinations Signal+Background,
203   // returns kFALSE when all combinations are made
204   // May be useful to introduce some options like "One-to-One", "All-to-One" and "All-to-All" ?
205
206   //realizing "One-to-One" option...
207
208   if(!fInitialized)
209     Init(1) ;
210
211   Int_t inputs ;
212   Bool_t endNotReached = kTRUE ;
213
214   for(inputs = 0; (inputs < fNinputs) && endNotReached ;inputs++){
215     if(fIevent->At(inputs)+1 < fIeventMax->At(inputs))
216       fIevent->AddAt(fIevent->At(inputs)+1, inputs ) ;
217     else
218       endNotReached = kFALSE ;
219   }
220   return endNotReached ;
221
222 }
223
224 //____________________________________________________________________________
225 void AliPHOSDigitizer::Digitize(Option_t *option) { 
226
227   //Makes the digitization of the collected summable digits
228
229   if(!fInitialized)
230     Init(1) ;
231
232   //Collects all hits in the same active volume into digit
233   //if(option == "raw")    // add simulated data to row data -- to be implemented
234
235
236   fDigits->Clear() ;
237
238   AliPHOS * PHOS = (AliPHOS *) gAlice->GetDetector("PHOS") ;   
239   AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance( PHOS->GetGeometry()->GetName(), PHOS->GetGeometry()->GetTitle() );
240
241   //Making digits with noise, first EMC
242   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
243   
244   Int_t nCPV ;
245   Int_t nPPSD ;
246   Int_t absID ;
247   TString name      =  geom->GetName() ;
248   
249   if ( name == "IHEP" || name == "MIXT" )    
250     nCPV =nEMC + geom->GetNumberOfCPVPadsZ()*geom->GetNumberOfCPVPadsPhi()*
251       geom->GetNCPVModules()*geom->GetNumberOfCPVLayers() ;
252   else
253     nCPV = nEMC; 
254   
255   if ( name == "GPS2" || name == "MIXT" )    
256     nPPSD =nCPV+2*geom->GetNPPSDModules()*geom->GetNumberOfModulesPhi()*geom->GetNumberOfModulesZ()*
257       geom->GetNumberOfPadsPhi()*geom->GetNumberOfPadsZ() ;
258   else
259     nPPSD = nCPV; 
260   
261   for(absID = 1; absID <= nEMC; absID++){
262     Float_t noise = gRandom->Gaus(0., fPinNoise) ; 
263     new((*fDigits)[absID-1]) AliPHOSDigit( -1,absID,fSDigitizer->Digitize(noise) ) ;
264   }
265   
266   for(absID = nEMC+1; absID <= nCPV; absID++){
267     Float_t noise = gRandom->Gaus(0., fCPVNoise) ; 
268     new((*fDigits)[absID-1]) AliPHOSDigit( -1,absID,fSDigitizer->Digitize(noise) ) ;
269   }
270   
271   for(absID = nCPV+1; absID <= nPPSD; absID++){
272     Float_t noise = gRandom->Gaus(0., fPPSDNoise) ; 
273     new((*fDigits)[absID-1]) AliPHOSDigit( -1,absID,fSDigitizer->Digitize(noise) ) ;
274   }
275   
276
277   // Now look throught (unsorted) list of SDigits and add corresponding digits  
278   AliPHOSDigit *curSDigit ;
279   AliPHOSDigit *digit ;
280     
281   Int_t inputs;
282   for(inputs = 0; inputs< fNinputs ; inputs++){  //loop over (possible) merge sources
283     
284     TClonesArray * sdigits= (TClonesArray *)fSDigits->At(inputs) ;
285     Int_t isdigit ;
286     Int_t nSDigits = sdigits->GetEntries() ;     
287     for(isdigit=0;isdigit< nSDigits; isdigit++){
288       curSDigit = (AliPHOSDigit *)sdigits->At(isdigit) ;
289       if(inputs)                                       //Shift primaries for non-background sdigits
290         curSDigit->ShiftPrimary(inputs) ;
291       digit = (AliPHOSDigit *)fDigits->At(curSDigit->GetId() - 1);
292       *digit = *digit + *curSDigit ;
293     }  
294   }
295
296
297
298   //remove digits below thresholds
299   for(absID = 0; absID < nEMC ; absID++)
300     if(fSDigitizer->Calibrate(((AliPHOSDigit*)fDigits->At(absID))->GetAmp()) < fEMCDigitThreshold)
301       fDigits->RemoveAt(absID) ;
302   for(absID = nEMC; absID < nCPV ; absID++)
303     if(fSDigitizer->Calibrate(((AliPHOSDigit*)fDigits->At(absID))->GetAmp()) < fCPVDigitThreshold)
304       fDigits->RemoveAt(absID) ;
305   for(absID = nCPV; absID < nPPSD ; absID++)
306     if(fSDigitizer->Calibrate(((AliPHOSDigit *)fDigits->At(absID))->GetAmp()) < fPPSDDigitThreshold)
307       fDigits->RemoveAt(absID) ;
308   
309   fDigits->Compress() ;  
310   
311   Int_t ndigits = fDigits->GetEntries() ;
312   fDigits->Expand(ndigits) ;
313
314
315   //Set indexes in list of digits
316   Int_t i ;
317   for (i = 0 ; i < ndigits ; i++) { 
318     AliPHOSDigit * digit = (AliPHOSDigit *) fDigits->At(i) ; 
319     digit->SetIndexInList(i) ;     
320   }
321 }
322 //____________________________________________________________________________
323 void AliPHOSDigitizer::WriteDigits(){
324
325   //Made TreeD in the output file if necessary and writes digiths there.
326
327   gAlice->GetEvent(fIevent->At(fOutFileNumber)) ;  // Suitable only for One-To-One mixing
328   gAlice->SetEvent(fIevent->At(fOutFileNumber)) ;  // for all-to-all will produce a lot of branches in TreeD
329
330   if(gAlice->TreeD()==0)
331     gAlice->MakeTree("D") ;
332   
333   //Make branches in TreeD for digits and Digitizer
334   char branchname[20];
335   sprintf(branchname,"PHOS");  
336
337   Int_t bufferSize = 16000 ;
338   char * filename = 0;
339   if(!fDigitsFile.IsNull())
340     filename = (char*) fDigitsFile.Data() ; //ievent ;
341   else
342     if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name
343       filename = new char[30] ;
344       //        sprintf(filename,"PHOS.Digits%d.root",ievent) ;
345       sprintf(filename,"PHOS.Digits.root") ;
346     }
347     else
348       filename = 0 ;
349   
350   //Link digits  
351   gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,&fDigits,bufferSize,filename);  
352   //Link Digitizer
353   AliPHOSDigitizer * d = this ;
354   Int_t splitlevel = 0 ;
355   sprintf(branchname,"AliPHOSDigitizer");   
356   gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,"AliPHOSDigitizer",&d, bufferSize, splitlevel,filename); 
357
358   gAlice->TreeD()->Fill() ;
359    
360   gAlice->TreeD()->Write(0,kOverwrite) ;  
361 }
362
363 //____________________________________________________________________________
364 void AliPHOSDigitizer::Exec(Option_t *option) { 
365   //manager
366
367   if(!fInitialized)    Init(1) ;
368   
369   while(Combinator()){  
370     
371     if(!ReadSDigits()) //read sdigits event(s) evaluated by Combinator() from file(s)
372       return ;    
373     
374     Digitize(option) ; //Add prepared SDigits to digits and add the noise
375     WriteDigits() ;
376   }
377
378 //   //Close all opened files
379 //   Int_t input ;
380 //   for(input = 0; input < fNinputs ; input ++){
381 //     TFile * file = (TFile*) gROOT->GetFile(((TObjString *) fHeaderFiles->At(input))->GetString() ) ;
382 //     file->Close() ;
383 //   }
384
385   
386 }
387
388 //__________________________________________________________________
389 Bool_t AliPHOSDigitizer::ReadSDigits(){
390 // Reads summable digits from the opened files for the particular set of events given by fIevent
391
392   if(!fInitialized)    Init(1) ;
393
394   Int_t inputs ;
395   for(inputs = fNinputs-1; inputs >= 0; inputs --){
396
397     Int_t event = fIevent->At(inputs) ;
398
399     TFile * file = (TFile*) gROOT->GetFile(((TObjString *) fHeaderFiles->At(inputs))->GetString() ) ;
400     file->cd() ;
401
402     // Get SDigits Tree header from file
403     char treeName[20]; 
404     sprintf(treeName,"TreeS%d",event);
405     TTree * treeS = (TTree*)file->Get(treeName);
406    
407     if(treeS==0){
408       cout << "Error at AliPHOSDigitizer: no "<<treeName << "   in file " << file->GetName() << endl ;
409       cout << "Do nothing " << endl ;
410       return kFALSE ;
411     }
412
413     TBranch * sdigitsBranch = 0;
414     TBranch * sdigitizerBranch = 0;
415
416     TObjArray * branches = treeS->GetListOfBranches() ;
417     Int_t ibranch;
418     Bool_t phosNotFound = kTRUE ;
419     Bool_t sdigitizerNotFound = kTRUE ;
420   
421     for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
422
423       if(phosNotFound){
424         sdigitsBranch=(TBranch *) branches->At(ibranch) ;
425         if( ((TObjString*)fSDigitsFiles->At(inputs))->GetString().CompareTo(sdigitsBranch->GetFileName())==0 ){
426           if( strcmp(sdigitsBranch->GetName(),"PHOS") == 0) {
427             phosNotFound = kFALSE ;
428           }
429         }
430       }
431
432       if(sdigitizerNotFound){
433         sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
434         if( ((TObjString*)fSDigitsFiles->At(inputs))->GetString().CompareTo(sdigitizerBranch->GetFileName()) == 0){
435           if( strcmp(sdigitizerBranch->GetName(),"AliPHOSSDigitizer") == 0) {
436             sdigitizerNotFound = kFALSE ;
437           }
438         }
439       }
440       
441     }
442
443     if(sdigitizerNotFound || phosNotFound){
444       cout << "Can't find Branch with sdigits or SDigitizer in the file " ;
445       if( ((TObjString*)fSDigitsFiles->At(inputs))->GetString().IsNull() )
446         cout << file->GetName() << endl ;       
447       else
448         cout << ((TObjString*)fSDigitsFiles->At(inputs))->GetString().Data() << endl ;
449       cout << "Do nothing" <<endl  ;
450       return kFALSE ;
451     }
452
453
454     TClonesArray * sdigits = (TClonesArray*) fSDigits->At(inputs) ;  
455     sdigitsBranch->SetAddress(&sdigits) ;
456     
457     AliPHOSSDigitizer *sDigitizer = new AliPHOSSDigitizer();
458     sdigitizerBranch->SetAddress(&sDigitizer) ;
459     treeS->GetEvent(0) ;
460
461     if(fSDigitizer == 0)
462       fSDigitizer = sDigitizer ;
463     else
464       if(!((*fSDigitizer)==(*sDigitizer)) ){
465         cout << "ERROR: you are using sdigits made with different SDigitizers" << endl ;
466         fSDigitizer->Print("") ;
467         sDigitizer->Print("") ;
468         cout << "Do Nothing " << endl ;
469         return kFALSE ;
470       }
471     
472   }
473   
474
475   return kTRUE ;
476
477 }
478 //__________________________________________________________________
479 void AliPHOSDigitizer::MixWith(char* HeaderFile,Int_t isOutFile = 1, char* SDigitsFile =0){
480 //
481
482   if(!fInitialized)
483     if(isOutFile)
484       Init(0) ;     //Do not read gAlice from Background file
485     else
486       Init(1) ;     //read gAlice from background file
487
488   if(HeaderFile == 0){
489     cout << "Specify at least header file to merge"<< endl ;
490     return ;
491   }
492   
493   Int_t inputs ;
494   for(inputs = 0; inputs < fNinputs ; inputs++){
495     if(strcmp(((TObjString *)fHeaderFiles->At(inputs))->GetString(),HeaderFile) == 0 ){
496       if(SDigitsFile == 0){ 
497         if(((TObjString*)fSDigitsFiles->At(inputs))->GetString().CompareTo("")  == 0){
498           cout << "Entry already exists, do not add" << endl ;
499           return ;
500         }
501       }
502       else
503         if(((TObjString*)fSDigitsFiles->At(inputs))->GetString().CompareTo(SDigitsFile) == 0){
504         cout << "Entry already exists, do not add" << endl ;
505         return;
506       }
507     }   
508   }  
509   
510   fHeaderFiles->Expand(fNinputs+1) ;
511   new((*fHeaderFiles)[fNinputs]) TObjString(HeaderFile) ;
512
513   
514   TFile * file ;
515   if(isOutFile)
516     file = new TFile(((TObjString *) fHeaderFiles->At(fNinputs))->GetString(),"update") ;  
517   else
518     file = new TFile(((TObjString *) fHeaderFiles->At(fNinputs))->GetString()) ;  
519
520   file->cd() ;
521
522   fSDigitsFiles->Expand(fNinputs+1) ;
523   new((*fSDigitsFiles)[fNinputs]) TObjString(SDigitsFile) ;
524
525   fSDigits->Expand(fNinputs+1) ;
526   new((*fSDigits)[fNinputs]) TClonesArray("AliPHOSDigit",1000) ;
527
528   fIevent->Set(fNinputs+1) ;
529   fIevent->AddAt(-1, fNinputs) ;
530
531   fIeventMax->Set(fNinputs+1) ;  
532
533   if(isOutFile){
534     gAlice = (AliRun*) file->Get("gAlice") ;
535     fIeventMax->AddAt((Int_t) gAlice->TreeE()->GetEntries(), fNinputs );
536     fOutFileNumber = fNinputs ;
537   }
538   else{
539     TTree * te = (TTree *) file->Get("TE") ;
540     fIeventMax->AddAt((Int_t) te->GetEntries(), fNinputs );
541   }
542
543   fNinputs++ ;
544
545 }
546 //__________________________________________________________________
547 void AliPHOSDigitizer::Print(Option_t* option){
548
549   if(!fInitialized)    Init(1) ;
550
551   cout << "------------------- "<< GetName() << " -------------" << endl ;
552   cout << "Digitizing sDigits from file(s): " <<endl ;
553   Int_t input ;
554   for(input = 0; input < fNinputs ; input++) {
555     cout << "          " << ((TObjString *) fHeaderFiles->At(input))->GetString() << 
556       "   Branch: " << ((TObjString *) fSDigitsFiles->At(input))->GetString() << endl ;
557   }
558   cout << endl ;
559   cout << "Writing digits to " << ((TObjString *) fHeaderFiles->At(fOutFileNumber))->GetString() << endl ;
560
561   cout << endl ;
562   cout << "With following parameters: " << endl ;
563   cout << "     Electronics noise in EMC (fPinNoise) = " << fPinNoise << endl ;
564   cout << "  Threshold  in EMC  (fEMCDigitThreshold) = " << fEMCDigitThreshold  << endl ; ;
565   cout << "                 Noise in CPV (fCPVNoise) = " << fCPVNoise << endl ; 
566   cout << "    Threshold in CPV (fCPVDigitThreshold) = " << fCPVDigitThreshold << endl ; 
567   cout << "               Noise in PPSD (fPPSDNoise) = " << fPPSDNoise << endl ;
568   cout << "  Threshold in PPSD (fPPSDDigitThreshold) = " << fPPSDDigitThreshold << endl ;
569   cout << "---------------------------------------------------" << endl ;
570  
571
572
573 }