Corrected for Bugs and added getter methods
[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"
814ad4bf 67#include "AliEMCALGetter.h"
61e0abb5 68
69ClassImp(AliEMCALSDigitizer)
70
71
72//____________________________________________________________________________
73 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
74{
75 // ctor
76 fA = 0;
e1f60236 77 fB = 0. ;
78 fPrimThreshold = 0. ;
61e0abb5 79 fNevents = 0 ;
80 fSDigits = 0 ;
81 fHits = 0 ;
e1f60236 82 fLayerRatio = 0. ;
61e0abb5 83 fIsInitialized = kFALSE ;
84
85}
86
87//____________________________________________________________________________
814ad4bf 88AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask(sDigitsTitle, headerFile)
61e0abb5 89{
90 // ctor
91 fA = 0;
92 fB = 10000000.;
ffa6d63b 93 fPrimThreshold = 0.01 ;
814ad4bf 94 fNevents = 0 ;
e1f60236 95 fLayerRatio = 5./6. ;
814ad4bf 96 Init();
61e0abb5 97}
98
99//____________________________________________________________________________
100AliEMCALSDigitizer::~AliEMCALSDigitizer()
101{
102 // dtor
103 if(fSDigits)
104 delete fSDigits ;
105 if(fHits)
106 delete fHits ;
107}
108//____________________________________________________________________________
109void AliEMCALSDigitizer::Init(){
61e0abb5 110
814ad4bf 111 // Initialization: open root-file, allocate arrays for hits and sdigits,
112 // attach task SDigitizer to the list of PHOS tasks
113 //
114 // Initialization can not be done in the default constructor
115 //============================================================= YS
116 // The initialisation is now done by the getter
61e0abb5 117
814ad4bf 118if( strcmp(GetTitle(), "") == 0 )
119 SetTitle("galice.root") ;
120
121 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), "update") ;
122 if ( gime == 0 ) {
123 cerr << "ERROR: AliEMCALSDigitizer::Init -> Could not obtain the Getter object !" << endl ;
124 return ;
125 }
126
127 gime->PostSDigits( GetName(), GetTitle() ) ;
61e0abb5 128
814ad4bf 129 TString sdname(GetName() );
130 sdname.Append(":") ;
131 sdname.Append(GetTitle() ) ;
132 SetName(sdname) ;
133 gime->PostSDigitizer(this) ;
134
135 // create a folder on the white board //YSAlice/WhiteBoard/SDigits/EMCAL/headerFile/sdigitsTitle
136
137 }
61e0abb5 138//____________________________________________________________________________
139void AliEMCALSDigitizer::Exec(Option_t *option) {
61e0abb5 140
814ad4bf 141// Collects all hits in the same active volume into digit
142
143 if( strcmp(GetName(), "") == 0 )
144 Init() ;
61e0abb5 145
814ad4bf 146 if (strstr(option, "print") ) {
147 Print("") ;
148 return ;
149 }
61e0abb5 150
814ad4bf 151 if(strstr(option,"tim"))
152 gBenchmark->Start("EMCALSDigitizer");
153 //Check, if this branch already exits
154 gAlice->GetEvent(0) ;
155 if(gAlice->TreeS() ) {
156 TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
157 TIter next(lob) ;
158 TBranch * branch = 0 ;
159 Bool_t emcalfound = kFALSE, sdigitizerfound = kFALSE ;
61e0abb5 160
814ad4bf 161 while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
162 if ( (strcmp(branch->GetName(), "EMCAL")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
163 emcalfound = kTRUE ;
164
165 else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
166 sdigitizerfound = kTRUE ;
61e0abb5 167 }
168
814ad4bf 169 if ( emcalfound || sdigitizerfound ) {
170 cerr << "WARNING: AliEMCALSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName()
171 << " already exits" << endl ;
172 return ;
173 }
174 }
175 TString sdname(GetName()) ;
176 sdname.Remove(sdname.Index(GetTitle())-1) ;
177 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
814ad4bf 178 Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
179
180 Int_t ievent ;
181 for(ievent = 0; ievent < nevents; ievent++){
182 gime->Event(ievent,"H") ;
183 const TClonesArray * fHits = gime->Hits() ;
184 TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
185 sdigits->Clear();
61e0abb5 186 Int_t nSdigits = 0 ;
814ad4bf 187
188 //Collects all hits in the same active volume into digit
189
190
61e0abb5 191
192
193 //Now made SDigits from hits, for EMCAL it is the same, so just copy
814ad4bf 194 // Int_t itrack ;
195 // for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
196 //gime->Track(itrack);
61e0abb5 197 //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
814ad4bf 198
199
200
61e0abb5 201
202 Int_t i;
e1f60236 203 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { // loop over all hits (hit = deposited energy/layer/entering particle)
204 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(fHits->At(i)) ;
205 AliEMCALDigit * curSDigit = 0 ;
206 AliEMCALDigit * sdigit = 0 ;
ffa6d63b 207 Bool_t newsdigit = kTRUE;
e1f60236 208
209 // Assign primary number only if deposited energy is significant
210
211 Float_t preShowerFactor ;
212 if (hit->IsInPreShower() )
213 preShowerFactor = fLayerRatio ;
214 else
215 preShowerFactor = 1 ;
216
217 if( hit->GetEnergy() > fPrimThreshold)
814ad4bf 218 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
e1f60236 219 hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
220 Digitize( hit->GetEnergy() * preShowerFactor ) ,
814ad4bf 221 hit->GetTime()) ;
61e0abb5 222 else
814ad4bf 223 curSDigit = new AliEMCALDigit( -1 ,
224 -1 ,
e1f60236 225 Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
226 Digitize( hit->GetEnergy() * preShowerFactor ) ,
227 hit->GetTime() ) ;
228 Int_t check = 0 ;
229 for(check= 0; check < nSdigits ; check++) {
814ad4bf 230 sdigit = (AliEMCALDigit *)sdigits->At(check);
e1f60236 231 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same tower or the same preshower ?
232 *sdigit = *sdigit + *curSDigit;
233
234 newsdigit = kFALSE;
235 }
814ad4bf 236 }
e1f60236 237 if (newsdigit) {
238 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
ffa6d63b 239 nSdigits++ ;
e1f60236 240 }
241 } // loop over all hits (hit = deposited energy/layer/entering particle)
814ad4bf 242 // } // loop over tracks
e1f60236 243
814ad4bf 244 sdigits->Sort() ;
61e0abb5 245
814ad4bf 246 nSdigits = sdigits->GetEntriesFast() ;
247 sdigits->Expand(nSdigits) ;
248
249 // Int_t i ;
e1f60236 250 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
251
252 Int_t lastPreShowerIndex = nSdigits - 1 ;
253 Int_t firstPreShowerIndex = -1 ;
254 Int_t index ;
255 AliEMCALDigit * sdigit = 0 ;
256 for ( index = 0; index < nSdigits ; index++) {
257 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ) ;
258 if (sdigit->IsInPreShower() ){
259 firstPreShowerIndex = index ;
260 break ;
261 }
262 }
263
264 AliEMCALDigit * preshower ;
265 AliEMCALDigit * tower ;
266 Int_t lastIndex = lastPreShowerIndex +1 ;
267
268
269 for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) {
270 preshower = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) );
271 Bool_t towerFound = kFALSE ;
272 Int_t jndex ;
273 for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
274 tower = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) );
275 if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {
276 *tower = *tower + *preshower ; // and add preshower to tower
277 towerFound = kTRUE ;
278 }
279 }
280 if ( !towerFound ) {
281
282 new((*sdigits)[lastIndex]) AliEMCALDigit(*preshower);
283 AliEMCALDigit * temp = dynamic_cast<AliEMCALDigit *>(sdigits->At(lastIndex)) ;
284 temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ;
285 lastIndex++ ;
286 }
287 }
288
289 sdigits->Sort() ;
290 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
291 sdigit->SetIndexInList(i) ;
61e0abb5 292 }
814ad4bf 293
294 if(gAlice->TreeS() == 0)
295 gAlice->MakeTree("S") ;
296
297
298
299
300 //Make (if necessary) branches
301 char * file =0;
302 if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
303 file = new char[strlen(gAlice->GetBaseFile())+20] ;
304 sprintf(file,"%s/EMCAL.SDigits.root",gAlice->GetBaseFile()) ;
61e0abb5 305 }
61e0abb5 306
814ad4bf 307 TDirectory *cwd = gDirectory;
308
309 //First list of sdigits
310 Int_t bufferSize = 32000 ;
311 TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
312 sdigitsBranch->SetTitle(sdname);
313 cout << " AliEMCALSDigitizer::Exec sdname " << sdname << endl ;
314
315 if (file) {
316 sdigitsBranch->SetFile(file);
317 TIter next( sdigitsBranch->GetListOfBranches());
318 TBranch * subbr;
319 while ((subbr=(TBranch*)next())) {
320 subbr->SetFile(file);
321 }
322 cwd->cd();
323 }
324
325 //second - SDigitizer
326 Int_t splitlevel = 0 ;
327 AliEMCALSDigitizer * sd = this ;
328 TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
329 &sd,bufferSize,splitlevel);
330 sdigitizerBranch->SetTitle(sdname);
331 if (file) {
332 sdigitizerBranch->SetFile(file);
333 TIter next( sdigitizerBranch->GetListOfBranches());
334 TBranch * subbr ;
335 while ((subbr=(TBranch*)next())) {
336 subbr->SetFile(file);
337 }
338 cwd->cd();
339 delete file;
340 }
341
342 sdigitsBranch->Fill() ;
343 sdigitizerBranch->Fill() ;
344 gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
345
346 if(strstr(option,"deb"))
347 PrintSDigits(option) ;
348
61e0abb5 349 }
350
351 if(strstr(option,"tim")){
352 gBenchmark->Stop("EMCALSDigitizer");
353 cout << "AliEMCALSDigitizer:" << endl ;
354 cout << " took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing "
355 << gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
356 cout << endl ;
357 }
358
359
360}
361//__________________________________________________________________
362void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
814ad4bf 363
364 // Setting title to branch SDigits
365
366 TString stitle(title) ;
367
368 // check if branch with title already exists
369 TBranch * sdigitsBranch =
370 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
371 TBranch * sdigitizerBranch =
372 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
373 const char * sdigitsTitle = sdigitsBranch ->GetTitle() ;
374 const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
375 if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
376 cerr << "ERROR: AliEMCALSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
377 return ;
378 }
379
380 cout << "AliEMCALSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
381
382 SetName(title) ;
383
384 // Post to the WhiteBoard
385 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
386 gime->PostSDigits( title, GetTitle()) ;
387
388
61e0abb5 389}
390//__________________________________________________________________
391void AliEMCALSDigitizer::Print(Option_t* option)const{
392 cout << "------------------- "<< GetName() << " -------------" << endl ;
814ad4bf 393 cout << " Writing SDigitis to branch with title " << GetName() << endl ;
61e0abb5 394 cout << " with digitization parameters A = " << fA << endl ;
395 cout << " B = " << fB << endl ;
396 cout << " Threshold for Primary assignment= " << fPrimThreshold << endl ;
397 cout << "---------------------------------------------------"<<endl ;
398
399}
400//__________________________________________________________________
401Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
402 if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
403 return kTRUE ;
404 else
405 return kFALSE ;
406}
407//__________________________________________________________________
408void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
409 //Prints list of digits produced at the current pass of AliEMCALDigitizer
410
814ad4bf 411 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
412 TString sdname(GetName()) ;
413 sdname.Remove(sdname.Index(GetTitle())-1) ;
414 const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
61e0abb5 415
814ad4bf 416 cout << "AliEMCALSDigitiser: event " << gAlice->GetEvNumber() << endl ;
417 cout << " Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
418 cout << endl ;
419 if(strstr(option,"all")||strstr(option,"EMC")){
420
61e0abb5 421 //loop over digits
422 AliEMCALDigit * digit;
e1f60236 423 cout << "SDigit Id " << " Amplitude " << " Time " << " Index " << " Nprim " << " Primaries list " << endl;
61e0abb5 424 Int_t index ;
814ad4bf 425 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
426 digit = (AliEMCALDigit * ) sdigits->At(index) ;
427 cout << setw(6) << digit->GetId() << " " << setw(10) << digit->GetAmp() << " " << digit->GetTime()
428 << setw(6) << digit->GetIndexInList() << " "
429 << setw(5) << digit->GetNprimary() <<" ";
61e0abb5 430
431 Int_t iprimary;
432 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
814ad4bf 433 cout << " " << digit->GetPrimary(iprimary+1) << " ";
61e0abb5 434 cout << endl;
435 }
814ad4bf 436 cout <<endl;
61e0abb5 437 }
438}
814ad4bf 439//________________________________________________________________________
440Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t s, Bool_t preshower)
441{
442 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
443 const AliEMCALGeometry * geom = gime->EMCALGeometry();
444 Int_t a = (s/geom->GetNPhi())%geom->GetNZ()+1; // Phi Tower Index
445 Int_t b = (s-1)%(geom->GetNPhi())+1; //Eta Tower Index
446 Int_t t = -10;
447
448if (a > 0 && b > 0)
449{
450 t = a*b + preshower*geom->GetNPhi()*geom->GetNZ();
451 return t;
452}
453else
454{cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Phi number = "
455 << b << "Eta number = " << a << endl ;
456return t;}
457}