Change to the directory of the file when the getter is created (fFile->cd())
[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"
106fc2fa 56#include "TGeometry.h"
61e0abb5 57#include "TBenchmark.h"
58// --- Standard library ---
59#include <iomanip.h>
60
61// --- AliRoot header files ---
62#include "AliRun.h"
106fc2fa 63#include "AliHeader.h"
61e0abb5 64#include "AliEMCALDigit.h"
65#include "AliEMCALHit.h"
66#include "AliEMCALSDigitizer.h"
67#include "AliEMCALGeometry.h"
68#include "AliEMCALv1.h"
814ad4bf 69#include "AliEMCALGetter.h"
61e0abb5 70
71ClassImp(AliEMCALSDigitizer)
72
73
74//____________________________________________________________________________
75 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
76{
77 // ctor
89e103bd 78 fA = fB = fNevents = 0 ;
79 fTowerPrimThreshold = fPreShowerPrimThreshold = fPhotonElectronFactor = 0. ;
80 fHits = fSDigits = fSDigits = 0 ;
81
61e0abb5 82 fIsInitialized = kFALSE ;
83
89e103bd 84
61e0abb5 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.;
89e103bd 93 fTowerPrimThreshold = 0.01 ;
94 fPreShowerPrimThreshold = 0.0001 ;
814ad4bf 95 fNevents = 0 ;
89e103bd 96 fPhotonElectronFactor = 5000. ; // photoelectrons per GeV
814ad4bf 97 Init();
61e0abb5 98}
99
100//____________________________________________________________________________
101AliEMCALSDigitizer::~AliEMCALSDigitizer()
102{
103 // dtor
104 if(fSDigits)
105 delete fSDigits ;
106 if(fHits)
107 delete fHits ;
108}
109//____________________________________________________________________________
110void AliEMCALSDigitizer::Init(){
61e0abb5 111
814ad4bf 112 // Initialization: open root-file, allocate arrays for hits and sdigits,
113 // attach task SDigitizer to the list of PHOS tasks
114 //
115 // Initialization can not be done in the default constructor
116 //============================================================= YS
117 // The initialisation is now done by the getter
61e0abb5 118
814ad4bf 119if( strcmp(GetTitle(), "") == 0 )
120 SetTitle("galice.root") ;
121
122 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), "update") ;
123 if ( gime == 0 ) {
124 cerr << "ERROR: AliEMCALSDigitizer::Init -> Could not obtain the Getter object !" << endl ;
125 return ;
126 }
127
128 gime->PostSDigits( GetName(), GetTitle() ) ;
61e0abb5 129
814ad4bf 130 TString sdname(GetName() );
131 sdname.Append(":") ;
132 sdname.Append(GetTitle() ) ;
97f3344e 133 SetName(sdname.Data()) ;
814ad4bf 134 gime->PostSDigitizer(this) ;
135
814ad4bf 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");
97f3344e 153
154
814ad4bf 155 //Check, if this branch already exits
156 gAlice->GetEvent(0) ;
97f3344e 157
158 TString sdname(GetName()) ;
159 sdname.Remove(sdname.Index(GetTitle())-1) ;
160
814ad4bf 161 if(gAlice->TreeS() ) {
162 TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
163 TIter next(lob) ;
164 TBranch * branch = 0 ;
165 Bool_t emcalfound = kFALSE, sdigitizerfound = kFALSE ;
61e0abb5 166
97f3344e 167 while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
814ad4bf 168 if ( (strcmp(branch->GetName(), "EMCAL")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
169 emcalfound = kTRUE ;
170
97f3344e 171 else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && (strcmp(branch->GetTitle(), sdname)==0) )
814ad4bf 172 sdigitizerfound = kTRUE ;
61e0abb5 173 }
174
814ad4bf 175 if ( emcalfound || sdigitizerfound ) {
176 cerr << "WARNING: AliEMCALSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName()
177 << " already exits" << endl ;
178 return ;
179 }
180 }
97f3344e 181
814ad4bf 182 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
814ad4bf 183 Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
184
185 Int_t ievent ;
186 for(ievent = 0; ievent < nevents; ievent++){
187 gime->Event(ievent,"H") ;
188 const TClonesArray * fHits = gime->Hits() ;
189 TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
190 sdigits->Clear();
61e0abb5 191 Int_t nSdigits = 0 ;
814ad4bf 192
193 //Collects all hits in the same active volume into digit
194
195
61e0abb5 196
197
198 //Now made SDigits from hits, for EMCAL it is the same, so just copy
814ad4bf 199 // Int_t itrack ;
200 // for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
201 //gime->Track(itrack);
61e0abb5 202 //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
89e103bd 203
61e0abb5 204 Int_t i;
e1f60236 205 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { // loop over all hits (hit = deposited energy/layer/entering particle)
206 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(fHits->At(i)) ;
207 AliEMCALDigit * curSDigit = 0 ;
208 AliEMCALDigit * sdigit = 0 ;
ffa6d63b 209 Bool_t newsdigit = kTRUE;
e1f60236 210
89e103bd 211
212
e1f60236 213 // Assign primary number only if deposited energy is significant
214
89e103bd 215 if( (!hit->IsInPreShower() && hit->GetEnergy() > fTowerPrimThreshold) ||
216 (hit->IsInPreShower() && hit->GetEnergy() > fPreShowerPrimThreshold)) {
217 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
e1f60236 218 hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
89e103bd 219 Digitize(hit->GetEnergy()),
814ad4bf 220 hit->GetTime()) ;
89e103bd 221 } else {
814ad4bf 222 curSDigit = new AliEMCALDigit( -1 ,
223 -1 ,
e1f60236 224 Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
89e103bd 225 Digitize(hit->GetEnergy()),
e1f60236 226 hit->GetTime() ) ;
89e103bd 227 }
e1f60236 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() ;
106fc2fa 247 if (nSdigits > 0) {
248 sdigits->Expand(nSdigits) ;
249
250 // Int_t i ;
251 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
252
253 Int_t lastPreShowerIndex = nSdigits - 1 ;
254 if (!(dynamic_cast<AliEMCALDigit *>(sdigits->At(lastPreShowerIndex))->IsInPreShower()))
255 lastPreShowerIndex = -2;
256 Int_t firstPreShowerIndex = 100000 ;
257 Int_t index ;
258 AliEMCALDigit * sdigit = 0 ;
259 for ( index = 0; index < nSdigits ; index++) {
260 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ) ;
261 if (sdigit->IsInPreShower() ){
262 firstPreShowerIndex = index ;
263 break ;
264 }
265 }
266
267 AliEMCALDigit * preshower ;
268 AliEMCALDigit * tower ;
269 Int_t lastIndex = lastPreShowerIndex +1 ;
270
271
272 for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) {
273 preshower = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) );
274 Bool_t towerFound = kFALSE ;
275 Int_t jndex ;
276 for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
277 tower = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) );
278 if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {
279 Float_t towerEnergy = static_cast<Float_t>(tower->GetAmp()) ;
280 Float_t preshoEnergy = static_cast<Float_t>(preshower->GetAmp()) ;
281 towerEnergy +=preshoEnergy ;
282 *tower = *tower + *preshower ; // and add preshower multiplied by layer ratio to tower
283 tower->SetAmp(static_cast<Int_t>(TMath::Ceil(towerEnergy))) ;
284 towerFound = kTRUE ;
285 }
286 }
287 if ( !towerFound ) {
288
289 new((*sdigits)[lastIndex]) AliEMCALDigit(*preshower);
290 AliEMCALDigit * temp = dynamic_cast<AliEMCALDigit *>(sdigits->At(lastIndex)) ;
291 temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ;
292 lastIndex++ ;
293 }
294 }
295
296 sdigits->Sort() ;
297 Int_t NPrimarymax = -1 ;
298 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
299 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
300 sdigit->SetIndexInList(i) ;
301 }
302
303 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
304 if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
305 NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
e1f60236 306 }
307 }
106fc2fa 308 if(gAlice->TreeS() == 0)
97f3344e 309 gAlice->MakeTree("S") ;
814ad4bf 310
814ad4bf 311 //First list of sdigits
312 Int_t bufferSize = 32000 ;
313 TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
314 sdigitsBranch->SetTitle(sdname);
814ad4bf 315
814ad4bf 316 //second - SDigitizer
317 Int_t splitlevel = 0 ;
318 AliEMCALSDigitizer * sd = this ;
319 TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
320 &sd,bufferSize,splitlevel);
321 sdigitizerBranch->SetTitle(sdname);
814ad4bf 322
ed99616a 323 gAlice->TreeS()->Fill() ;
106fc2fa 324 gAlice->TreeS()->AutoSave() ;
814ad4bf 325
326 if(strstr(option,"deb"))
327 PrintSDigits(option) ;
328
61e0abb5 329 }
330
331 if(strstr(option,"tim")){
332 gBenchmark->Stop("EMCALSDigitizer");
333 cout << "AliEMCALSDigitizer:" << endl ;
334 cout << " took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing "
335 << gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
336 cout << endl ;
337 }
338
339
340}
341//__________________________________________________________________
342void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
814ad4bf 343
344 // Setting title to branch SDigits
345
346 TString stitle(title) ;
347
348 // check if branch with title already exists
349 TBranch * sdigitsBranch =
350 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
351 TBranch * sdigitizerBranch =
352 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
353 const char * sdigitsTitle = sdigitsBranch ->GetTitle() ;
354 const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
355 if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
356 cerr << "ERROR: AliEMCALSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
357 return ;
358 }
359
360 cout << "AliEMCALSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
361
362 SetName(title) ;
363
364 // Post to the WhiteBoard
365 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
366 gime->PostSDigits( title, GetTitle()) ;
367
368
61e0abb5 369}
106fc2fa 370
371//__________________________________________________________________
372void AliEMCALSDigitizer::SetSplitFile(const TString splitFileName) const
373{
374 // Diverts the SDigits in a file separate from the hits file
375
376 TDirectory * cwd = gDirectory ;
377 TFile * splitFile = gAlice->InitTreeFile("S",splitFileName.Data());
378 splitFile->cd() ;
379 gAlice->Write();
380
381 TTree *treeE = gAlice->TreeE();
382 if (!treeE) {
383 cerr<<"No TreeE found "<<endl;
384 abort() ;
385 }
386
387 // copy TreeE
388 AliHeader *header = new AliHeader();
389 treeE->SetBranchAddress("Header", &header);
390 treeE->SetBranchStatus("*",1);
391 TTree *treeENew = treeE->CloneTree();
392 treeENew->Write();
393
394 // copy AliceGeom
395 TGeometry *AliceGeom = static_cast<TGeometry*>(cwd->Get("AliceGeom"));
396 if (!AliceGeom) {
397 cerr<<"AliceGeom was not found in the input file "<<endl;
398 abort() ;
399 }
400 AliceGeom->Write();
401 cwd->cd() ;
402 gAlice->MakeTree("S",splitFile);
403 cout << "INFO: AliEMCALSDigitizer::SetSPlitMode -> SDigits will be stored in " << splitFileName.Data() << endl ;
404}
405
61e0abb5 406//__________________________________________________________________
407void AliEMCALSDigitizer::Print(Option_t* option)const{
408 cout << "------------------- "<< GetName() << " -------------" << endl ;
814ad4bf 409 cout << " Writing SDigitis to branch with title " << GetName() << endl ;
89e103bd 410 cout << " with digitization parameters A = " << fA << endl ;
411 cout << " B = " << fB << endl ;
412 cout << " Threshold for Primary assignment in Tower = " << fTowerPrimThreshold << endl ;
413 cout << " Threshold for Primary assignment in PreShower = " << fPreShowerPrimThreshold << endl ;
61e0abb5 414 cout << "---------------------------------------------------"<<endl ;
415
416}
417//__________________________________________________________________
418Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
89e103bd 419 if( (fA==sd.fA)&&(fB==sd.fB)&&
420 (fTowerPrimThreshold==sd.fTowerPrimThreshold) &&
421 (fPreShowerPrimThreshold==sd.fPreShowerPrimThreshold))
61e0abb5 422 return kTRUE ;
423 else
424 return kFALSE ;
425}
426//__________________________________________________________________
427void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
428 //Prints list of digits produced at the current pass of AliEMCALDigitizer
429
814ad4bf 430 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
431 TString sdname(GetName()) ;
432 sdname.Remove(sdname.Index(GetTitle())-1) ;
433 const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
61e0abb5 434
814ad4bf 435 cout << "AliEMCALSDigitiser: event " << gAlice->GetEvNumber() << endl ;
436 cout << " Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
437 cout << endl ;
438 if(strstr(option,"all")||strstr(option,"EMC")){
439
61e0abb5 440 //loop over digits
441 AliEMCALDigit * digit;
e1f60236 442 cout << "SDigit Id " << " Amplitude " << " Time " << " Index " << " Nprim " << " Primaries list " << endl;
61e0abb5 443 Int_t index ;
814ad4bf 444 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
445 digit = (AliEMCALDigit * ) sdigits->At(index) ;
446 cout << setw(6) << digit->GetId() << " " << setw(10) << digit->GetAmp() << " " << digit->GetTime()
447 << setw(6) << digit->GetIndexInList() << " "
448 << setw(5) << digit->GetNprimary() <<" ";
61e0abb5 449
450 Int_t iprimary;
451 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
814ad4bf 452 cout << " " << digit->GetPrimary(iprimary+1) << " ";
61e0abb5 453 cout << endl;
454 }
814ad4bf 455 cout <<endl;
61e0abb5 456 }
457}
814ad4bf 458//________________________________________________________________________
229f77c4 459Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower){
460 // Method to Transform from Hit Id to Digit Id
461 // This function should be one to one
814ad4bf 462 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
463 const AliEMCALGeometry * geom = gime->EMCALGeometry();
f063936c 464 Int_t ieta = ((ihit-1)/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index
229f77c4 465 Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index
466 Int_t it = -10;
467 Int_t ipre = 0;
814ad4bf 468
229f77c4 469 if (preshower)ipre = 1;
470 if (iphi > 0 && ieta >= 0){
471 it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ();
472 return it;
473 }else{
474 cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Eta number = "
475 << ieta << "Phi number = " << iphi << endl ;
476 return it;
477 } // end if iphi>0 && ieta>0
814ad4bf 478}
229f77c4 479//_______________________________________________________________________________________
89e103bd 480// void AliEMCALSDigitizer::TestTowerID(void)
481// {
482// Int_t j;
229f77c4 483
89e103bd 484// Bool_t preshower = kFALSE;
485// for (j = 0 ; j < 10 ; j++){ // loop over hit id
486// Int_t i;
487// for (i = 0 ; i <= 2 ; i++){ // loop over
488// Int_t k = i*96*144+j*144+1;
489// cout << " Hit Index = " << k << " " << j*10 << " TOWERID = " << Layer2TowerID(k, preshower) << endl ;
490// }
491// }
492// }