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