Adapted to PHOSGetter
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
CommitLineData
ffa6d63b 1/*************************************************************************
61e0abb5 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 makes SDigits out of Hits
20// A Summable Digits is the sum of all hits originating
ffa6d63b 21// from one in one tower of the EMCAL
61e0abb5 22// A threshold for assignment of the primary to SDigit is applied
23// SDigits are written to TreeS, branch "EMCAL"
24// AliEMCALSDigitizer with all current parameters is written
25// to TreeS branch "AliEMCALSDigitizer".
26// Both branches have the same title. If necessary one can produce
27// another set of SDigits with different parameters. Two versions
28// can be distunguished using titles of the branches.
29// User case:
30// root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root")
31// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
32// root [1] s->ExecuteTask()
33// // Makes SDigitis for all events stored in galice.root
34// root [2] s->SetPedestalParameter(0.001)
35// // One can change parameters of digitization
36// root [3] s->SetSDigitsBranch("Redestal 0.001")
37// // and write them into the new branch
38// root [4] s->ExeciteTask("deb all tim")
39// // available parameters:
40// deb - print # of produced SDigitis
41// deb all - print # and list of produced SDigits
42// tim - print benchmarking information
43//
ffa6d63b 44//*-- Author : Sahal Yacoob (LBL)
45// based on : AliPHOSSDigitzer
61e0abb5 46//////////////////////////////////////////////////////////////////////////////
47
48
49// --- ROOT system ---
50#include "TFile.h"
51#include "TTask.h"
52#include "TTree.h"
53#include "TSystem.h"
54#include "TROOT.h"
55#include "TFolder.h"
56#include "TBenchmark.h"
57// --- Standard library ---
58#include <iomanip.h>
59
60// --- AliRoot header files ---
61#include "AliRun.h"
62#include "AliEMCALDigit.h"
63#include "AliEMCALHit.h"
64#include "AliEMCALSDigitizer.h"
65#include "AliEMCALGeometry.h"
66#include "AliEMCALv1.h"
67
68ClassImp(AliEMCALSDigitizer)
69
70
71//____________________________________________________________________________
72 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
73{
74 // ctor
75 fA = 0;
76 fB = 10000000. ;
ffa6d63b 77 fPrimThreshold = 0.01 ;
61e0abb5 78 fNevents = 0 ;
79 fSDigits = 0 ;
80 fHits = 0 ;
81 fIsInitialized = kFALSE ;
82
83}
84
85//____________________________________________________________________________
86AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliEMCALSDigitizer","")
87{
88 // ctor
89 fA = 0;
90 fB = 10000000.;
ffa6d63b 91 fPrimThreshold = 0.01 ;
61e0abb5 92 fNevents = 0 ;
93 fSDigitsTitle = sDigitsTitle ;
94 fHeadersFile = headerFile ;
d75bea67 95 fSDigits = new TClonesArray("AliEMCALDigit",30000);
61e0abb5 96 fHits = new TClonesArray("AliEMCALHit",1000);
97
98 TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
99
100 //File was not opened yet
101 if(file == 0){
102 if(fHeadersFile.Contains("rfio"))
103 file = TFile::Open(fHeadersFile,"update") ;
104 else
105 file = new TFile(fHeadersFile.Data(),"update") ;
106 gAlice = (AliRun *) file->Get("gAlice") ;
107 }
108
109 //add Task to //root/Tasks folder
110 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
111 roottasks->Add(this) ;
112
113 fIsInitialized = kTRUE ;
114}
115
116//____________________________________________________________________________
117AliEMCALSDigitizer::~AliEMCALSDigitizer()
118{
119 // dtor
120 if(fSDigits)
121 delete fSDigits ;
122 if(fHits)
123 delete fHits ;
124}
125//____________________________________________________________________________
126void AliEMCALSDigitizer::Init(){
127 //Initialization can not be done in the default constructor
128
129 if(!fIsInitialized){
130
131 if(fHeadersFile.IsNull())
132 fHeadersFile="galice.root" ;
133
134 TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
135
136 //if file was not opened yet, read gAlice
137 if(file == 0){
138 file = new TFile(fHeadersFile.Data(),"update") ;
139 gAlice = (AliRun *) file->Get("gAlice") ;
140 }
141
142 fHits = new TClonesArray("AliEMCALHit",1000);
d75bea67 143 fSDigits = new TClonesArray("AliEMCALDigit",30000);
61e0abb5 144
145 // add Task to //root/Tasks folder
146 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
147 roottasks->Add(this) ;
148
149 fIsInitialized = kTRUE ;
150 }
151}
152//____________________________________________________________________________
153void AliEMCALSDigitizer::Exec(Option_t *option) {
154 //Collects all hits in the same active volume into digit
155
156 if(!fIsInitialized)
157 Init() ;
158
159 if(strstr(option,"tim"))
160 gBenchmark->Start("EMCALSDigitizer");
161
162 fNevents = (Int_t) gAlice->TreeE()->GetEntries() ;
163
164 Int_t ievent ;
165 for(ievent = 0; ievent < fNevents; ievent++){
166 gAlice->GetEvent(ievent) ;
167 gAlice->SetEvent(ievent) ;
168
169 if(gAlice->TreeH()==0){
170 cout << "AliEMCALSDigitizer: There is no Hit Tree" << endl;
171 return ;
172 }
173
174 //set address of the hits
175 TBranch * branch = gAlice->TreeH()->GetBranch("EMCAL");
176 if (branch)
177 branch->SetAddress(&fHits);
178 else{
179 cout << "ERROR in AliEMCALSDigitizer: "<< endl ;
180 cout << " no branch EMCAL in TreeH"<< endl ;
181 cout << " do nothing " << endl ;
182 return ;
183 }
184
185 fSDigits->Clear();
186 Int_t nSdigits = 0 ;
187
188
189 //Now made SDigits from hits, for EMCAL it is the same, so just copy
190 Int_t itrack ;
191 for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
192
193 //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
194 branch->GetEntry(itrack,0);
ffa6d63b 195 AliEMCALv0 * EMCAL = (AliEMCALv0 *) gAlice->GetDetector("EMCAL") ;
61e0abb5 196 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance( EMCAL->GetGeometry()->GetName(), EMCAL->GetGeometry()->GetTitle() );
197
198 Int_t i;
199 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
200 AliEMCALHit * hit = (AliEMCALHit*)fHits->At(i) ;
ffa6d63b 201 AliEMCALDigit * curSDigit = 0 ;
202 AliEMCALDigit * sdigit ;
203 Bool_t newsdigit = kTRUE;
204// Assign primary number only if contribution is significant
61e0abb5 205 if( hit->GetEnergy() > fPrimThreshold)
d75bea67 206 curSDigit = new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1 ) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
61e0abb5 207 else
d75bea67 208 curSDigit = new AliEMCALDigit( -1 , -1 ,(((hit->GetId()/geom->GetNPhi())%geom->GetNZ() + 1 ) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
209 //cout << "Hit ID = " <<hit->GetId() << endl ;
210 //cout << "ID for detector = " << curSDigit->GetId() << endl ;
211 // cout << hit->GetEnergy() << " - hit energy - Digit Energy - " << curSDigit->GetAmp() << endl;
ffa6d63b 212 for(Int_t check= 0; check < nSdigits ; check++) {
61e0abb5 213 sdigit = (AliEMCALDigit *)fSDigits->At(check);
ffa6d63b 214 if( sdigit->GetId() == curSDigit->GetId())
d75bea67 215 {// cout << "SDigit - Get Amp " << sdigit->GetAmp() << endl ;
ffa6d63b 216 *sdigit = *sdigit + *curSDigit ;
217 newsdigit = kFALSE;
d75bea67 218 // cout << " and after addition " << sdigit->GetAmp() << endl ;
61e0abb5 219 }
ffa6d63b 220 }
221 if (newsdigit)
222 { new((*fSDigits)[nSdigits]) AliEMCALDigit(*curSDigit);
223 nSdigits++ ;
d75bea67 224 //cout << "Detector nsdigits = " << nSdigits << endl ;
225 }
ffa6d63b 226 newsdigit = kTRUE;
227
61e0abb5 228
229 if( hit->GetEnergy() > fPrimThreshold)
d75bea67 230 curSDigit = new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), ((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ() + 1) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
61e0abb5 231 else
d75bea67 232 curSDigit = new AliEMCALDigit( -1 , -1 ,((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * ((hit->GetId()-1)%(geom->GetNPhi())+1)), Digitize( hit->GetEnergy() ) ) ;
61e0abb5 233
234 if((hit->GetId()/geom->GetNPhi()) < (2*geom->GetNZ()))
d75bea67 235 { //cout << "ID for Preshower = " << curSDigit->GetId() << endl ;
61e0abb5 236 for(Int_t check= 0; check < nSdigits; check++) {
237 sdigit = (AliEMCALDigit *)fSDigits->At(check);
ffa6d63b 238 if( sdigit->GetId() == curSDigit->GetId())
61e0abb5 239 {
ffa6d63b 240 *sdigit = *sdigit + *curSDigit ;
241 newsdigit = kFALSE ;
242 }
61e0abb5 243 }
ffa6d63b 244 if (newsdigit)
245 { new((*fSDigits)[nSdigits]) AliEMCALDigit(*curSDigit);
246 nSdigits++ ;
d75bea67 247 //cout << "Preshower nsdigits = " << nSdigits << endl ;
248 }
ffa6d63b 249 newsdigit=kTRUE;
61e0abb5 250 }
251 }
252 } // loop over tracks
61e0abb5 253 fSDigits->Sort() ;
254
255 nSdigits = fSDigits->GetEntriesFast() ;
256 fSDigits->Expand(nSdigits) ;
257
258 Int_t i ;
259 for (i = 0 ; i < nSdigits ; i++) {
260 AliEMCALDigit * digit = (AliEMCALDigit *) fSDigits->At(i) ;
261 digit->SetIndexInList(i) ;
262 }
263
264 if(gAlice->TreeS() == 0)
265 gAlice->MakeTree("S") ;
266
267 //check, if this branch already exits?
268 TBranch * sdigitsBranch = 0;
269 TBranch * sdigitizerBranch = 0;
270
271 TObjArray * branches = gAlice->TreeS()->GetListOfBranches() ;
272 Int_t ibranch;
273 Bool_t emcalNotFound = kTRUE ;
274 Bool_t sdigitizerNotFound = kTRUE ;
275
276 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
277
278 if(emcalNotFound){
279 sdigitsBranch=(TBranch *) branches->At(ibranch) ;
280 if( (strcmp("EMCAL",sdigitsBranch->GetName())==0 ) &&
281 (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) )
282 emcalNotFound = kFALSE ;
283 }
284 if(sdigitizerNotFound){
285 sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
286 if( (strcmp(sdigitizerBranch->GetName(),"AliEMCALSDigitizer") == 0)&&
287 (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) )
288 sdigitizerNotFound = kFALSE ;
289 }
290 }
291
292 if(!(sdigitizerNotFound && emcalNotFound)){
293 cout << "AliEMCALSdigitizer error:" << endl ;
294 cout << "Can not overwrite existing branches: do not write" << endl ;
295 return ;
296 }
297
298 //Make (if necessary) branches
299 char * file =0;
300 if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
301 file = new char[strlen(gAlice->GetBaseFile())+20] ;
302 sprintf(file,"%s/EMCAL.SDigits.root",gAlice->GetBaseFile()) ;
303 }
304
305 TDirectory *cwd = gDirectory;
306
307 //First list of sdigits
308 Int_t bufferSize = 32000 ;
309 sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&fSDigits,bufferSize);
310 sdigitsBranch->SetTitle(fSDigitsTitle.Data());
311 if (file) {
312 sdigitsBranch->SetFile(file);
313 TIter next( sdigitsBranch->GetListOfBranches());
314 TBranch * subbr;
315 while ((subbr=(TBranch*)next())) {
316 subbr->SetFile(file);
317 }
318 cwd->cd();
319 }
320
321 //second - SDigitizer
322 Int_t splitlevel = 0 ;
323 AliEMCALSDigitizer * sd = this ;
324 sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
325 &sd,bufferSize,splitlevel);
326 sdigitizerBranch->SetTitle(fSDigitsTitle.Data());
327 if (file) {
328 sdigitizerBranch->SetFile(file);
329 TIter next( sdigitizerBranch->GetListOfBranches());
330 TBranch * subbr ;
331 while ((subbr=(TBranch*)next())) {
332 subbr->SetFile(file);
333 }
334 cwd->cd();
335 delete file;
336 }
337
338 sdigitsBranch->Fill() ;
339 sdigitizerBranch->Fill() ;
340 gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
341
342 if(strstr(option,"deb"))
343 PrintSDigits(option) ;
344
345 }
346
347 if(strstr(option,"tim")){
348 gBenchmark->Stop("EMCALSDigitizer");
349 cout << "AliEMCALSDigitizer:" << endl ;
350 cout << " took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing "
351 << gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
352 cout << endl ;
353 }
354
355
356}
357//__________________________________________________________________
358void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
359 //Seting title to branch SDigits
360 if(!fSDigitsTitle.IsNull())
361 cout << "AliEMCALSdigitizer: changing SDigits file from " <<fSDigitsTitle.Data() << " to " << title << endl ;
362 fSDigitsTitle=title ;
363}
364//__________________________________________________________________
365void AliEMCALSDigitizer::Print(Option_t* option)const{
366 cout << "------------------- "<< GetName() << " -------------" << endl ;
367 cout << " Writing SDigitis to branch with title " << fSDigitsTitle.Data() << endl ;
368 cout << " with digitization parameters A = " << fA << endl ;
369 cout << " B = " << fB << endl ;
370 cout << " Threshold for Primary assignment= " << fPrimThreshold << endl ;
371 cout << "---------------------------------------------------"<<endl ;
372
373}
374//__________________________________________________________________
375Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
376 if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
377 return kTRUE ;
378 else
379 return kFALSE ;
380}
381//__________________________________________________________________
382void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
383 //Prints list of digits produced at the current pass of AliEMCALDigitizer
384
385 cout << "AliEMCALSDigitizer: " << endl ;
386 cout << " Number of entries in SDigits list " << fSDigits->GetEntriesFast() << endl ;
387 cout << endl ;
388
389 if(strstr(option,"all")){// print all digits
390
391 //loop over digits
392 AliEMCALDigit * digit;
393 cout << "SDigit Id " << " Amplitude " << " Index " << " Nprim " << " Primaries list " << endl;
394 Int_t index ;
395 for (index = 0 ; index < fSDigits->GetEntries() ; index++) {
396 digit = (AliEMCALDigit * ) fSDigits->At(index) ;
397 cout << setw(8) << digit->GetId() << " " << setw(3) << digit->GetAmp() << " "
398 << setw(6) << digit->GetIndexInList() << " "
399 << setw(5) << digit->GetNprimary() <<" ";
400
401 Int_t iprimary;
402 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
403 cout << setw(5) << digit->GetPrimary(iprimary+1) << " ";
404 cout << endl;
405 }
406
407 }
408}