added the HCAL section and a method to get the local position
[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
05a92d59 46// Modif:
47// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
48// of new IO (à la PHOS)
61e0abb5 49//////////////////////////////////////////////////////////////////////////////
50
61e0abb5 51// --- ROOT system ---
52#include "TFile.h"
53#include "TTask.h"
54#include "TTree.h"
55#include "TSystem.h"
56#include "TROOT.h"
57#include "TFolder.h"
58#include "TBenchmark.h"
05a92d59 59#include "TGeometry.h"
60
61e0abb5 61// --- Standard library ---
61e0abb5 62
63// --- AliRoot header files ---
64#include "AliRun.h"
106fc2fa 65#include "AliHeader.h"
61e0abb5 66#include "AliEMCALDigit.h"
61e0abb5 67#include "AliEMCALGeometry.h"
814ad4bf 68#include "AliEMCALGetter.h"
05a92d59 69#include "AliEMCALHit.h"
70#include "AliEMCALSDigitizer.h"
61e0abb5 71
72ClassImp(AliEMCALSDigitizer)
61e0abb5 73
74//____________________________________________________________________________
75 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
76{
77 // ctor
839828a6 78 InitParameters() ;
92f521a9 79 fDefaultInit = kTRUE ;
61e0abb5 80}
81
82//____________________________________________________________________________
05a92d59 83AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle, const Bool_t toSplit):
84TTask(sDigitsTitle, headerFile)
61e0abb5 85{
86 // ctor
839828a6 87 InitParameters() ;
05a92d59 88 fToSplit = toSplit ;
814ad4bf 89 Init();
92f521a9 90 fDefaultInit = kFALSE ;
61e0abb5 91}
92
93//____________________________________________________________________________
94AliEMCALSDigitizer::~AliEMCALSDigitizer()
95{
96 // dtor
839828a6 97 fSplitFile = 0 ;
61e0abb5 98}
36657114 99
61e0abb5 100//____________________________________________________________________________
101void AliEMCALSDigitizer::Init(){
814ad4bf 102 // Initialization: open root-file, allocate arrays for hits and sdigits,
36657114 103 // attach task SDigitizer to the list of EMCAL tasks
814ad4bf 104 //
105 // Initialization can not be done in the default constructor
106 //============================================================= YS
107 // The initialisation is now done by the getter
61e0abb5 108
36657114 109 if( strcmp(GetTitle(), "") == 0 )
814ad4bf 110 SetTitle("galice.root") ;
afa0c708 111
05a92d59 112 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), fToSplit) ;
814ad4bf 113 if ( gime == 0 ) {
9859bfc0 114 Error("Init", "Could not obtain the Getter object !" ) ;
814ad4bf 115 return ;
116 }
117
118 gime->PostSDigits( GetName(), GetTitle() ) ;
05a92d59 119
120 fSplitFile = 0 ;
121 if(fToSplit){
122 // construct the name of the file as /path/EMCAL.SDigits.root
123 // First - extract full path if necessary
124 TString sDigitsFileName(GetTitle()) ;
125 Ssiz_t islash = sDigitsFileName.Last('/') ;
126 if(islash<sDigitsFileName.Length())
127 sDigitsFileName.Remove(islash+1,sDigitsFileName.Length()) ;
128 else
129 sDigitsFileName="" ;
173558f2 130
05a92d59 131 // Next - append the file name
132 sDigitsFileName+="EMCAL.SDigits." ;
133 if((strcmp(GetName(),"Default")!=0)&&(strcmp(GetName(),"")!=0)){
134 sDigitsFileName+=GetName() ;
135 sDigitsFileName+="." ;
136 }
137 sDigitsFileName+="root" ;
173558f2 138
05a92d59 139 // Finally - check if the file already opened or open the file
140 fSplitFile = static_cast<TFile*>(gROOT->GetFile(sDigitsFileName.Data()));
141 if(!fSplitFile)
142 fSplitFile = TFile::Open(sDigitsFileName.Data(),"update") ;
143 }
61e0abb5 144
814ad4bf 145 TString sdname(GetName() );
146 sdname.Append(":") ;
147 sdname.Append(GetTitle() ) ;
05a92d59 148 SetName(sdname) ;
814ad4bf 149 gime->PostSDigitizer(this) ;
839828a6 150}
151
152//____________________________________________________________________________
afa0c708 153void AliEMCALSDigitizer::InitParameters()
154{
05a92d59 155 fA = 0;
156 fB = 10000000.;
157 fTowerPrimThreshold = 0.01 ;
839828a6 158 fPreShowerPrimThreshold = 0.0001 ;
05a92d59 159 fPhotonElectronFactor = 5000. ; // photoelectrons per GeV
160 fSplitFile = 0 ;
161 fToSplit = kFALSE ;
839828a6 162}
163
61e0abb5 164//____________________________________________________________________________
afa0c708 165void AliEMCALSDigitizer::Exec(Option_t *option)
166{
167 // Collects all hits in the same active volume into digit
173558f2 168
814ad4bf 169 if( strcmp(GetName(), "") == 0 )
170 Init() ;
61e0abb5 171
814ad4bf 172 if (strstr(option, "print") ) {
173 Print("") ;
174 return ;
175 }
61e0abb5 176
814ad4bf 177 if(strstr(option,"tim"))
178 gBenchmark->Start("EMCALSDigitizer");
97f3344e 179
814ad4bf 180 //Check, if this branch already exits
05a92d59 181 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
182 if(gime->BranchExists("SDigits") )
183 return;
184
97f3344e 185 TString sdname(GetName()) ;
186 sdname.Remove(sdname.Index(GetTitle())-1) ;
187
05a92d59 188 Int_t nevents = gime->MaxEvent() ;
814ad4bf 189 Int_t ievent ;
afa0c708 190
191 for(ievent = 0; ievent < nevents; ievent++){
192 gime->Event(ievent,"H") ;
193
194 const TClonesArray * hits = gime->Hits() ;
195
814ad4bf 196 TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
197 sdigits->Clear();
61e0abb5 198 Int_t nSdigits = 0 ;
814ad4bf 199
05a92d59 200 //Collects all hits in the same active volume into digit
61e0abb5 201
05a92d59 202 //Now make SDigits from hits, for EMCAL it is the same, so just copy
203 Int_t nPrim = static_cast<Int_t>((gAlice->TreeH())->GetEntries()) ;
204 // Attention nPrim is the number of primaries tracked by Geant
205 // and this number could be different to the number of Primaries in TreeK;
206 Int_t iprim ;
207 for ( iprim = 0 ; iprim < nPrim ; iprim++ ) {
208 //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
209 gime->Track(iprim) ;
210 Int_t i;
211 for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
212 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
213 AliEMCALDigit * curSDigit = 0 ;
214 AliEMCALDigit * sdigit = 0 ;
215 Bool_t newsdigit = kTRUE;
216 // Assign primary number only if deposited energy is significant
e1f60236 217
05a92d59 218 if( (!hit->IsInPreShower() && hit->GetEnergy() > fTowerPrimThreshold) ||
219 (hit->IsInPreShower() && hit->GetEnergy() > fPreShowerPrimThreshold))
220 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
221 hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
222 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
223 else
224 curSDigit = new AliEMCALDigit( -1 ,
225 -1 ,
226 Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
227 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
228 Int_t check = 0 ;
229 for(check= 0; check < nSdigits ; check++) {
230 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
231 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same tower or the same preshower ?
232 *sdigit = *sdigit + *curSDigit;
e1f60236 233 newsdigit = kFALSE;
05a92d59 234 }
e1f60236 235 }
05a92d59 236 if (newsdigit) {
237 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
238 nSdigits++ ;
239 }
afa0c708 240 delete curSDigit ;
05a92d59 241 } // loop over all hits (hit = deposited energy/layer/entering particle)
242 } // loop over iprim
e1f60236 243
814ad4bf 244 sdigits->Sort() ;
61e0abb5 245
814ad4bf 246 nSdigits = sdigits->GetEntriesFast() ;
05a92d59 247 fSDigitsInRun += nSdigits ;
248 sdigits->Expand(nSdigits) ;
106fc2fa 249
05a92d59 250 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
251
afa0c708 252 if (nSdigits != 0 ) {
253 Int_t lastPreShowerIndex = nSdigits - 1 ;
254
255
256 if (!(dynamic_cast<AliEMCALDigit *>(sdigits->At(lastPreShowerIndex))->IsInPreShower()))
257
258 lastPreShowerIndex = -2;
259
260 Int_t firstPreShowerIndex = 100000 ;
261 Int_t index ;
262 AliEMCALDigit * sdigit = 0 ;
263
264 for ( index = 0; index < nSdigits ; index++) {
265 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ) ;
266 if (sdigit->IsInPreShower() ){
267 firstPreShowerIndex = index ;
268 break ;
269 }
106fc2fa 270 }
271
afa0c708 272 AliEMCALDigit * preshower ;
273 AliEMCALDigit * tower ;
274 Int_t lastIndex = lastPreShowerIndex +1 ;
275
276 for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) {
277 preshower = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) );
278 Bool_t towerFound = kFALSE ;
279 Int_t jndex ;
280 for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
281 tower = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) );
282 if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {
283 Float_t towerEnergy = static_cast<Float_t>(tower->GetAmp()) ;
284 Float_t preshoEnergy = static_cast<Float_t>(preshower->GetAmp()) ;
285 towerEnergy +=preshoEnergy ;
286 *tower = *tower + *preshower ; // and add preshower multiplied by layer ratio to tower
287 tower->SetAmp(static_cast<Int_t>(TMath::Ceil(towerEnergy))) ;
288 towerFound = kTRUE ;
289 }
290 }
291 if ( !towerFound ) {
292 new((*sdigits)[lastIndex]) AliEMCALDigit(*preshower);
293 AliEMCALDigit * temp = dynamic_cast<AliEMCALDigit *>(sdigits->At(lastIndex)) ;
294 temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ;
295 lastIndex++ ;
106fc2fa 296 }
297 }
afa0c708 298 sdigits->Sort() ;
299 Int_t NPrimarymax = -1 ;
300 Int_t i ;
301 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
302 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
303 sdigit->SetIndexInList(i) ;
304 }
305
306 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
307 if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
308 NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
e1f60236 309 }
05a92d59 310 }
311
05a92d59 312 //Now write SDigits
313
314 if(gAlice->TreeS() == 0 || (fSplitFile)) //<--- To be checked: we should not create TreeS if it is already here
315 gAlice->MakeTree("S",fSplitFile);
afa0c708 316
05a92d59 317 if(fSplitFile)
318 fSplitFile->cd() ;
319
814ad4bf 320 //First list of sdigits
321 Int_t bufferSize = 32000 ;
322 TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
323 sdigitsBranch->SetTitle(sdname);
814ad4bf 324
05a92d59 325 //NEXT - SDigitizer
814ad4bf 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);
814ad4bf 331
36657114 332 sdigitsBranch->Fill() ;
333 sdigitizerBranch->Fill() ;
106fc2fa 334 gAlice->TreeS()->AutoSave() ;
05a92d59 335
814ad4bf 336 if(strstr(option,"deb"))
337 PrintSDigits(option) ;
338
61e0abb5 339 }
05a92d59 340
61e0abb5 341 if(strstr(option,"tim")){
9859bfc0 342 gBenchmark->Stop("EMCALSDigitizer");
343 Info("Exec", "took %f seconds for SDigitizing %f seconds per event",
344 gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ;
05a92d59 345 }
61e0abb5 346}
05a92d59 347
61e0abb5 348//__________________________________________________________________
349void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
814ad4bf 350
351 // Setting title to branch SDigits
352
353 TString stitle(title) ;
354
355 // check if branch with title already exists
356 TBranch * sdigitsBranch =
357 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
358 TBranch * sdigitizerBranch =
359 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
360 const char * sdigitsTitle = sdigitsBranch ->GetTitle() ;
361 const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
362 if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
9859bfc0 363 Error("SetSDigitsBranch", "Cannot overwrite existing branch with title %s", title) ;
814ad4bf 364 return ;
365 }
366
9859bfc0 367 Info("SetSDigitsBranch", "Changing SDigits file from %s to %s", GetName(), title) ;
814ad4bf 368
369 SetName(title) ;
370
371 // Post to the WhiteBoard
372 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
373 gime->PostSDigits( title, GetTitle()) ;
61e0abb5 374}
106fc2fa 375
106fc2fa 376
61e0abb5 377//__________________________________________________________________
05a92d59 378void AliEMCALSDigitizer::Print(Option_t* option)const
379{
380 // Prints parameters of SDigitizer
381
9859bfc0 382 TString message("\n") ;
383 message += "------------------- ";
384 message += GetName() ;
385 message += " -------------\n" ;
386 message += " Writing SDigitis to branch with title " ;
387 message += GetName() ;
388 message += "\n with digitization parameters A = " ;
389 message += fA ;
390 message += "\n B = " ;
391 message += fB ;
392 message += "\n Threshold for Primary assignment in Tower = " ;
393 message += fTowerPrimThreshold ;
394 message += "\n Threshold for Primary assignment in PreShower = " ;
395 message += fPreShowerPrimThreshold ;
396 message += "\n---------------------------------------------------" ;
61e0abb5 397
9859bfc0 398 Info("Print", message.Data() ) ;
61e0abb5 399}
05a92d59 400
61e0abb5 401//__________________________________________________________________
05a92d59 402Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
403{
404 // Equal operator.
405 // SDititizers are equal if their pedestal, slope and threshold are equal
406
89e103bd 407 if( (fA==sd.fA)&&(fB==sd.fB)&&
408 (fTowerPrimThreshold==sd.fTowerPrimThreshold) &&
409 (fPreShowerPrimThreshold==sd.fPreShowerPrimThreshold))
61e0abb5 410 return kTRUE ;
411 else
412 return kFALSE ;
413}
173558f2 414
61e0abb5 415//__________________________________________________________________
416void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
417 //Prints list of digits produced at the current pass of AliEMCALDigitizer
418
814ad4bf 419 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
420 TString sdname(GetName()) ;
421 sdname.Remove(sdname.Index(GetTitle())-1) ;
422 const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
61e0abb5 423
9859bfc0 424 TString message("\n") ;
425 message += "event " ;
426 message += gAlice->GetEvNumber() ;
427 message += "\n Number of entries in SDigits list " ;
428 message += sdigits->GetEntriesFast() ;
814ad4bf 429 if(strstr(option,"all")||strstr(option,"EMC")){
9859bfc0 430
61e0abb5 431 //loop over digits
432 AliEMCALDigit * digit;
11f9c5ff 433 message += "\n Id Amplitude Time Index Nprim: Primaries list \n" ;
61e0abb5 434 Int_t index ;
11f9c5ff 435 char * tempo = new char[8192];
814ad4bf 436 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
05a92d59 437 digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
11f9c5ff 438 sprintf(tempo, "\n%6d %8d %6.5e %4d %2d :",
439 digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
440 message += tempo ;
61e0abb5 441
442 Int_t iprimary;
9859bfc0 443 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
11f9c5ff 444 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
445 message += tempo ;
9859bfc0 446 }
61e0abb5 447 }
11f9c5ff 448 delete tempo ;
61e0abb5 449 }
9859bfc0 450 Info("PrintSDigits", message.Data() ) ;
61e0abb5 451}
05a92d59 452
814ad4bf 453//________________________________________________________________________
05a92d59 454const Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower)
455{
229f77c4 456 // Method to Transform from Hit Id to Digit Id
457 // This function should be one to one
814ad4bf 458 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
459 const AliEMCALGeometry * geom = gime->EMCALGeometry();
f063936c 460 Int_t ieta = ((ihit-1)/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index
229f77c4 461 Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index
462 Int_t it = -10;
463 Int_t ipre = 0;
814ad4bf 464
229f77c4 465 if (preshower)ipre = 1;
466 if (iphi > 0 && ieta >= 0){
467 it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ();
468 return it;
469 }else{
9859bfc0 470 Error("Layer2TowerID", "there is an error: Eta number = %f Phi number = %f", ieta, iphi) ;
229f77c4 471 return it;
472 } // end if iphi>0 && ieta>0
814ad4bf 473}
229f77c4 474//_______________________________________________________________________________________
89e103bd 475// void AliEMCALSDigitizer::TestTowerID(void)
476// {
477// Int_t j;
229f77c4 478
89e103bd 479// Bool_t preshower = kFALSE;
480// for (j = 0 ; j < 10 ; j++){ // loop over hit id
481// Int_t i;
482// for (i = 0 ; i <= 2 ; i++){ // loop over
483// Int_t k = i*96*144+j*144+1;
9859bfc0 484// Info("TestTowerID", " Hit Index = %d %d TOWERID = %d", k, j*10, Layer2TowerID(k, preshower) ) ;
89e103bd 485// }
486// }
487// }
05a92d59 488
489//____________________________________________________________________________
490void AliEMCALSDigitizer::UseHitsFrom(const char * filename)
491{
492 SetTitle(filename) ;
493 Init() ;
494}