Adaption to new fluka common blocks (E. Futo)
[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() ;
62741b41 79 fSplitFile = 0 ;
80 fToSplit = kFALSE ;
92f521a9 81 fDefaultInit = kTRUE ;
61e0abb5 82}
83
84//____________________________________________________________________________
05a92d59 85AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle, const Bool_t toSplit):
86TTask(sDigitsTitle, headerFile)
61e0abb5 87{
88 // ctor
05a92d59 89 fToSplit = toSplit ;
814ad4bf 90 Init();
62741b41 91 InitParameters() ;
92f521a9 92 fDefaultInit = kFALSE ;
61e0abb5 93}
94
95//____________________________________________________________________________
96AliEMCALSDigitizer::~AliEMCALSDigitizer()
97{
98 // dtor
839828a6 99 fSplitFile = 0 ;
61e0abb5 100}
36657114 101
61e0abb5 102//____________________________________________________________________________
103void AliEMCALSDigitizer::Init(){
814ad4bf 104 // Initialization: open root-file, allocate arrays for hits and sdigits,
36657114 105 // attach task SDigitizer to the list of EMCAL tasks
814ad4bf 106 //
107 // Initialization can not be done in the default constructor
108 //============================================================= YS
109 // The initialisation is now done by the getter
61e0abb5 110
36657114 111 if( strcmp(GetTitle(), "") == 0 )
814ad4bf 112 SetTitle("galice.root") ;
afa0c708 113
05a92d59 114 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), fToSplit) ;
814ad4bf 115 if ( gime == 0 ) {
9859bfc0 116 Error("Init", "Could not obtain the Getter object !" ) ;
814ad4bf 117 return ;
118 }
119
120 gime->PostSDigits( GetName(), GetTitle() ) ;
05a92d59 121
122 fSplitFile = 0 ;
62741b41 123
05a92d59 124 if(fToSplit){
125 // construct the name of the file as /path/EMCAL.SDigits.root
126 // First - extract full path if necessary
127 TString sDigitsFileName(GetTitle()) ;
128 Ssiz_t islash = sDigitsFileName.Last('/') ;
129 if(islash<sDigitsFileName.Length())
130 sDigitsFileName.Remove(islash+1,sDigitsFileName.Length()) ;
131 else
132 sDigitsFileName="" ;
173558f2 133
05a92d59 134 // Next - append the file name
135 sDigitsFileName+="EMCAL.SDigits." ;
136 if((strcmp(GetName(),"Default")!=0)&&(strcmp(GetName(),"")!=0)){
137 sDigitsFileName+=GetName() ;
138 sDigitsFileName+="." ;
139 }
140 sDigitsFileName+="root" ;
173558f2 141
05a92d59 142 // Finally - check if the file already opened or open the file
143 fSplitFile = static_cast<TFile*>(gROOT->GetFile(sDigitsFileName.Data()));
144 if(!fSplitFile)
145 fSplitFile = TFile::Open(sDigitsFileName.Data(),"update") ;
146 }
61e0abb5 147
814ad4bf 148 TString sdname(GetName() );
149 sdname.Append(":") ;
150 sdname.Append(GetTitle() ) ;
05a92d59 151 SetName(sdname) ;
814ad4bf 152 gime->PostSDigitizer(this) ;
62741b41 153
839828a6 154}
155
156//____________________________________________________________________________
afa0c708 157void AliEMCALSDigitizer::InitParameters()
158{
05a92d59 159 fA = 0;
160 fB = 10000000.;
62741b41 161
162 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
163 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
164 if (geom->GetSampling() == 0.) {
165 Error("InitParameters", "Sampling factor not set !") ;
166 abort() ;
167 }
168 else
169 Info("InitParameters", "Sampling factor set to %f\n", geom->GetSampling()) ;
170
171 // this threshold corresponds approximately to 100 MeV
172 fECPrimThreshold = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNECLayers() ;
173 fPREPrimThreshold = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNPRLayers() ;
174 fHCPrimThreshold = fECPrimThreshold/5. ; // 5 is totally arbitrary
175
839828a6 176}
177
61e0abb5 178//____________________________________________________________________________
afa0c708 179void AliEMCALSDigitizer::Exec(Option_t *option)
180{
62741b41 181 // Collects all hits in the section (PRE/ECAL/HCAL) of the same tower into digit
173558f2 182
814ad4bf 183 if( strcmp(GetName(), "") == 0 )
184 Init() ;
61e0abb5 185
814ad4bf 186 if (strstr(option, "print") ) {
187 Print("") ;
188 return ;
189 }
61e0abb5 190
814ad4bf 191 if(strstr(option,"tim"))
192 gBenchmark->Start("EMCALSDigitizer");
97f3344e 193
814ad4bf 194 //Check, if this branch already exits
05a92d59 195 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
196 if(gime->BranchExists("SDigits") )
197 return;
198
97f3344e 199 TString sdname(GetName()) ;
200 sdname.Remove(sdname.Index(GetTitle())-1) ;
201
05a92d59 202 Int_t nevents = gime->MaxEvent() ;
62741b41 203 Int_t ievent ;
204 for(ievent = 0; ievent < nevents; ievent++){
205 gime->Event(ievent,"H") ;
206 const TClonesArray * hits = gime->Hits() ;
814ad4bf 207 TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
208 sdigits->Clear();
61e0abb5 209 Int_t nSdigits = 0 ;
61e0abb5 210
05a92d59 211 //Now make SDigits from hits, for EMCAL it is the same, so just copy
212 Int_t nPrim = static_cast<Int_t>((gAlice->TreeH())->GetEntries()) ;
213 // Attention nPrim is the number of primaries tracked by Geant
214 // and this number could be different to the number of Primaries in TreeK;
215 Int_t iprim ;
62741b41 216 for ( iprim = 0 ; iprim < nPrim ; iprim++ ) {
217 //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
218 gime->Track(iprim) ;
219 Int_t i;
220 for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
221 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
222 AliEMCALDigit * curSDigit = 0 ;
223 AliEMCALDigit * sdigit = 0 ;
224 Bool_t newsdigit = kTRUE;
225
226 // Assign primary number only if deposited energy is significant
227
228 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
229
bda276cf 230 if (gDebug) {
62741b41 231 Info("Exec", "id = %d energy = %f thresholdPRE = %f thresholdEC = %f thresholdHC = %f \n",
232 hit->GetId(), hit->GetEnergy(), fPREPrimThreshold, fECPrimThreshold, fHCPrimThreshold) ;
bda276cf 233 Info("Exec", "where PRE/ECAL/HCAL %d, %d, %d", geom->IsInPRE(hit->GetId()), geom->IsInECAL(hit->GetId()), geom->IsInHCAL(hit->GetId()) ) ;
234 }
62741b41 235 if( geom->IsInPRE(hit->GetId()) )
236 if( hit->GetEnergy() > fPREPrimThreshold )
05a92d59 237 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
62741b41 238 hit->GetIparent(), hit->GetId(),
05a92d59 239 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
62741b41 240 else
05a92d59 241 curSDigit = new AliEMCALDigit( -1 ,
242 -1 ,
62741b41 243 hit->GetId(),
244 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
245 else if( geom->IsInECAL(hit->GetId()) )
246 if( hit->GetEnergy() > fECPrimThreshold )
247 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
248 hit->GetIparent(), hit->GetId(),
249 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
250 else
251 curSDigit = new AliEMCALDigit( -1 ,
252 -1 ,
253 hit->GetId(),
254 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
255 else if( geom->IsInHCAL(hit->GetId()) )
256 if( hit->GetEnergy() > fHCPrimThreshold )
257
258 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
259 hit->GetIparent(), hit->GetId(),
260 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
261 else
262 curSDigit = new AliEMCALDigit( -1 ,
263 -1 ,
264 hit->GetId(),
265 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
afa0c708 266
62741b41 267 Int_t check = 0 ;
268 for(check= 0; check < nSdigits ; check++) {
269 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
270 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL/HCAL/preshower tower ?
271 *sdigit = *sdigit + *curSDigit;
272 newsdigit = kFALSE;
afa0c708 273 }
afa0c708 274 }
62741b41 275 if (newsdigit) {
276 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
277 nSdigits++ ;
e1f60236 278 }
62741b41 279 delete curSDigit ;
280 } // loop over all hits (hit = deposited energy/layer/entering particle)
281 } // loop over iprim
282
283 sdigits->Sort() ;
284
285 nSdigits = sdigits->GetEntriesFast() ;
286 fSDigitsInRun += nSdigits ;
287 sdigits->Expand(nSdigits) ;
288
289 Int_t NPrimarymax = -1 ;
290 Int_t i ;
291 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
292 AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
293 sdigit->SetIndexInList(i) ;
294 }
295
296 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
297 if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
298 NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
299 }
300
301 // Now write SDigits
302
303 if(gAlice->TreeS() == 0 || (fSplitFile)) //<--- To be checked: we should not create TreeS if it is already here
304 gAlice->MakeTree("S",fSplitFile);
305
306 if(fSplitFile)
307 fSplitFile->cd() ;
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
314 //NEXT - SDigitizer
315 Int_t splitlevel = 0 ;
316 AliEMCALSDigitizer * sd = this ;
317 TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
318 &sd,bufferSize,splitlevel);
319 sdigitizerBranch->SetTitle(sdname);
320
321 sdigitsBranch->Fill() ;
322 sdigitizerBranch->Fill() ;
323
324 gAlice->TreeS()->AutoSave() ;
325
326 if(strstr(option,"deb"))
327 PrintSDigits(option) ;
61e0abb5 328 }
62741b41 329
61e0abb5 330 if(strstr(option,"tim")){
9859bfc0 331 gBenchmark->Stop("EMCALSDigitizer");
332 Info("Exec", "took %f seconds for SDigitizing %f seconds per event",
62741b41 333 gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ;
334 }
61e0abb5 335}
05a92d59 336
61e0abb5 337//__________________________________________________________________
338void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
814ad4bf 339
340 // Setting title to branch SDigits
341
342 TString stitle(title) ;
343
344 // check if branch with title already exists
345 TBranch * sdigitsBranch =
346 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
347 TBranch * sdigitizerBranch =
348 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
349 const char * sdigitsTitle = sdigitsBranch ->GetTitle() ;
350 const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
351 if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
9859bfc0 352 Error("SetSDigitsBranch", "Cannot overwrite existing branch with title %s", title) ;
814ad4bf 353 return ;
354 }
355
9859bfc0 356 Info("SetSDigitsBranch", "Changing SDigits file from %s to %s", GetName(), title) ;
814ad4bf 357
358 SetName(title) ;
359
360 // Post to the WhiteBoard
361 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
362 gime->PostSDigits( title, GetTitle()) ;
61e0abb5 363}
106fc2fa 364
106fc2fa 365
61e0abb5 366//__________________________________________________________________
05a92d59 367void AliEMCALSDigitizer::Print(Option_t* option)const
368{
369 // Prints parameters of SDigitizer
370
9859bfc0 371 TString message("\n") ;
372 message += "------------------- ";
373 message += GetName() ;
374 message += " -------------\n" ;
375 message += " Writing SDigitis to branch with title " ;
376 message += GetName() ;
377 message += "\n with digitization parameters A = " ;
378 message += fA ;
379 message += "\n B = " ;
380 message += fB ;
9859bfc0 381 message += "\n Threshold for Primary assignment in PreShower = " ;
62741b41 382 message += fPREPrimThreshold ;
383 message += "\n Threshold for Primary assignment in EC section= " ;
384 message += fECPrimThreshold ;
385 message += "\n Threshold for Primary assignment in HC section= " ;
386 message += fHCPrimThreshold ;
9859bfc0 387 message += "\n---------------------------------------------------" ;
61e0abb5 388
9859bfc0 389 Info("Print", message.Data() ) ;
61e0abb5 390}
05a92d59 391
61e0abb5 392//__________________________________________________________________
05a92d59 393Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
394{
395 // Equal operator.
396 // SDititizers are equal if their pedestal, slope and threshold are equal
397
89e103bd 398 if( (fA==sd.fA)&&(fB==sd.fB)&&
62741b41 399 (fECPrimThreshold==sd.fECPrimThreshold) &&
400 (fHCPrimThreshold==sd.fHCPrimThreshold) &&
401 (fPREPrimThreshold==sd.fPREPrimThreshold))
61e0abb5 402 return kTRUE ;
403 else
404 return kFALSE ;
405}
173558f2 406
61e0abb5 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
9859bfc0 416 TString message("\n") ;
417 message += "event " ;
418 message += gAlice->GetEvNumber() ;
419 message += "\n Number of entries in SDigits list " ;
420 message += sdigits->GetEntriesFast() ;
814ad4bf 421 if(strstr(option,"all")||strstr(option,"EMC")){
9859bfc0 422
61e0abb5 423 //loop over digits
424 AliEMCALDigit * digit;
11f9c5ff 425 message += "\n Id Amplitude Time Index Nprim: Primaries list \n" ;
61e0abb5 426 Int_t index ;
11f9c5ff 427 char * tempo = new char[8192];
814ad4bf 428 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
05a92d59 429 digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
11f9c5ff 430 sprintf(tempo, "\n%6d %8d %6.5e %4d %2d :",
431 digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
432 message += tempo ;
61e0abb5 433
434 Int_t iprimary;
9859bfc0 435 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
11f9c5ff 436 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
437 message += tempo ;
9859bfc0 438 }
61e0abb5 439 }
11f9c5ff 440 delete tempo ;
61e0abb5 441 }
9859bfc0 442 Info("PrintSDigits", message.Data() ) ;
61e0abb5 443}
05a92d59 444
229f77c4 445//_______________________________________________________________________________________
89e103bd 446// void AliEMCALSDigitizer::TestTowerID(void)
447// {
448// Int_t j;
229f77c4 449
89e103bd 450// Bool_t preshower = kFALSE;
451// for (j = 0 ; j < 10 ; j++){ // loop over hit id
452// Int_t i;
453// for (i = 0 ; i <= 2 ; i++){ // loop over
454// Int_t k = i*96*144+j*144+1;
9859bfc0 455// Info("TestTowerID", " Hit Index = %d %d TOWERID = %d", k, j*10, Layer2TowerID(k, preshower) ) ;
89e103bd 456// }
457// }
458// }
05a92d59 459
460//____________________________________________________________________________
461void AliEMCALSDigitizer::UseHitsFrom(const char * filename)
462{
463 SetTitle(filename) ;
464 Init() ;
465}