Most PWG1 AddTask macros used in the pilot QA train were modified to use mgr->GetComm...
[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
ab6a174f 23//
24// JLK 26-Jun-2008 Added explanation:
25// SDigits need to hold the energy sum of the hits, but AliEMCALDigit
26// can (should) only store amplitude. Therefore, the SDigit energy is
27// "digitized" before being stored and must be "calibrated" back to an
28// energy before SDigits are summed to form true Digits
29//
61e0abb5 30// SDigits are written to TreeS, branch "EMCAL"
31// AliEMCALSDigitizer with all current parameters is written
32// to TreeS branch "AliEMCALSDigitizer".
33// Both branches have the same title. If necessary one can produce
34// another set of SDigits with different parameters. Two versions
35// can be distunguished using titles of the branches.
36// User case:
37// root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root")
38// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
39// root [1] s->ExecuteTask()
40// // Makes SDigitis for all events stored in galice.root
41// root [2] s->SetPedestalParameter(0.001)
42// // One can change parameters of digitization
43// root [3] s->SetSDigitsBranch("Redestal 0.001")
44// // and write them into the new branch
45// root [4] s->ExeciteTask("deb all tim")
46// // available parameters:
47// deb - print # of produced SDigitis
48// deb all - print # and list of produced SDigits
49// tim - print benchmarking information
50//
ffa6d63b 51//*-- Author : Sahal Yacoob (LBL)
52// based on : AliPHOSSDigitzer
61e0abb5 53//////////////////////////////////////////////////////////////////////////////
54
61e0abb5 55// --- ROOT system ---
bda7cf56 56#include <TBenchmark.h>
bda7cf56 57#include <TBrowser.h>
d804d556 58//#include <Riostream.h>
bda7cf56 59#include <TMath.h>
a1e17193 60#include <TROOT.h>
05a92d59 61
61e0abb5 62// --- Standard library ---
d64c959b 63#include "stdlib.h"
61e0abb5 64
65// --- AliRoot header files ---
13e1db42 66#include "AliLog.h"
5dee926e 67#include "AliRunLoader.h"
68#include "AliStack.h"
61e0abb5 69#include "AliEMCALDigit.h"
5dee926e 70#include "AliEMCALLoader.h"
05a92d59 71#include "AliEMCALHit.h"
72#include "AliEMCALSDigitizer.h"
88cb7938 73#include "AliEMCALGeometry.h"
33a52e55 74#include "AliEMCALSimParam.h"
61e0abb5 75
76ClassImp(AliEMCALSDigitizer)
61e0abb5 77
78//____________________________________________________________________________
18a21c7c 79AliEMCALSDigitizer::AliEMCALSDigitizer()
80 : TTask("",""),
81 fA(0.),fB(0.),fECPrimThreshold(0.),
82 fDefaultInit(kTRUE),
83 fEventFolderName(0),
84 fInit(0),
85 fSDigitsInRun(0),
86 fFirstEvent(0),
87 fLastEvent(0),
7ea6391b 88 fSampling(0.)
61e0abb5 89{
90 // ctor
0a4cb131 91 InitParameters();
61e0abb5 92}
93
94//____________________________________________________________________________
e191bb57 95AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName,
18a21c7c 96 const char * eventFolderName)
97 : TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
98 fA(0.),fB(0.),fECPrimThreshold(0.),
99 fDefaultInit(kFALSE),
100 fEventFolderName(eventFolderName),
101 fInit(0),
102 fSDigitsInRun(0),
103 fFirstEvent(0),
104 fLastEvent(0),
7ea6391b 105 fSampling(0.)
61e0abb5 106{
107 // ctor
814ad4bf 108 Init();
62741b41 109 InitParameters() ;
7ea6391b 110
61e0abb5 111}
112
88cb7938 113
61e0abb5 114//____________________________________________________________________________
18a21c7c 115AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd)
116 : TTask(sd.GetName(),sd.GetTitle()),
117 fA(sd.fA),
118 fB(sd.fB),
119 fECPrimThreshold(sd.fECPrimThreshold),
120 fDefaultInit(sd.fDefaultInit),
121 fEventFolderName(sd.fEventFolderName),
122 fInit(sd.fInit),
123 fSDigitsInRun(sd.fSDigitsInRun),
124 fFirstEvent(sd.fFirstEvent),
125 fLastEvent(sd.fLastEvent),
7ea6391b 126 fSampling(sd.fSampling)
18a21c7c 127{
88cb7938 128 //cpy ctor
61e0abb5 129}
36657114 130
88cb7938 131
61e0abb5 132//____________________________________________________________________________
797e311f 133AliEMCALSDigitizer::~AliEMCALSDigitizer() {
14ce0a6e 134 //dtor
5dee926e 135 AliLoader *emcalLoader=0;
33c3c91a 136 if ((emcalLoader = AliRunLoader::Instance()->GetDetectorLoader("EMCAL")))
5dee926e 137 emcalLoader->CleanSDigitizer();
797e311f 138}
139
140//____________________________________________________________________________
61e0abb5 141void AliEMCALSDigitizer::Init(){
814ad4bf 142 // Initialization: open root-file, allocate arrays for hits and sdigits,
36657114 143 // attach task SDigitizer to the list of EMCAL tasks
814ad4bf 144 //
145 // Initialization can not be done in the default constructor
146 //============================================================= YS
147 // The initialisation is now done by the getter
61e0abb5 148
dc7da436 149 fInit = kTRUE;
afa0c708 150
33c3c91a 151 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
5dee926e 152
153 if ( emcalLoader == 0 ) {
154 Fatal("Init", "Could not obtain the AliEMCALLoader");
814ad4bf 155 return ;
156 }
157
5dee926e 158 emcalLoader->PostSDigitizer(this);
159 emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
88cb7938 160
839828a6 161}
162
163//____________________________________________________________________________
afa0c708 164void AliEMCALSDigitizer::InitParameters()
165{
14ce0a6e 166 //initialize parameters for sdigitization
167
5dee926e 168 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
62741b41 169 if (geom->GetSampling() == 0.) {
5082063d 170 Fatal("InitParameters", "Sampling factor not set !") ;
62741b41 171 }
ab6a174f 172
6569f329 173 // Get the parameters from the OCDB via the loader
174 AliRunLoader *rl = AliRunLoader::Instance();
175 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
176 AliEMCALSimParam * simParam = 0x0;
177 if(emcalLoader) simParam = emcalLoader->SimulationParameters();
178
179 if(!simParam){
180 simParam = AliEMCALSimParam::GetInstance();
181 AliWarning("Simulation Parameters not available in OCDB?");
182 }
183
ab6a174f 184 //
185 //JLK 26-Jun-2008 THIS SHOULD HAVE BEEN EXPLAINED AGES AGO:
186 //
187 //In order to be able to store SDigit Energy info into
188 //AliEMCALDigit, we need to convert it temporarily to an ADC amplitude
189 //and later when summing SDigits to form digits, convert it back to
190 //energy. These fA and fB parameters accomplish this through the
191 //Digitize() and Calibrate() methods
192 //
1963b290 193 // Initializes parameters
6569f329 194 fA = simParam->GetA(); //0;
195 fB = simParam->GetB(); //1.e+6; // Changed 24 Apr 2007. Dynamic range now 2 TeV
1963b290 196 fSampling = geom->GetSampling();
197
ff98b189 198 // threshold for deposit energy of hit
6569f329 199 fECPrimThreshold = simParam->GetECPrimaryThreshold();//0.05;// GeV // 22-may-07 was 0// 24-nov-04 - was 1.e-6;
ff98b189 200
b42d4197 201 AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
202 AliDebug(2,Form(" fInit %i\n", int(fInit)));
203 AliDebug(2,Form(" fFirstEvent %i\n", fFirstEvent));
204 AliDebug(2,Form(" fLastEvent %i\n", fLastEvent));
205 AliDebug(2,Form(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()));
206 AliDebug(2,Form(" with digitization parameters A = %f\n", fA));
207 AliDebug(2,Form(" B = %f\n", fB));
208 AliDebug(2,Form(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold));
209 AliDebug(2,Form(" Sampling = %f\n", fSampling));
210 AliDebug(2,Form("---------------------------------------------------\n"));
211
839828a6 212}
213
61e0abb5 214//____________________________________________________________________________
afa0c708 215void AliEMCALSDigitizer::Exec(Option_t *option)
216{
fdebddeb 217 // Collects all hit of the same tower into digits
1963b290 218 TString o(option); o.ToUpper();
814ad4bf 219 if (strstr(option, "print") ) {
b42d4197 220
221 AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
222 AliDebug(2,Form(" fInit %i\n", int(fInit)));
223 AliDebug(2,Form(" fFirstEvent %i\n", fFirstEvent));
224 AliDebug(2,Form(" fLastEvent %i\n", fLastEvent));
225 AliDebug(2,Form(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()));
226 AliDebug(2,Form(" with digitization parameters A = %f\n", fA));
227 AliDebug(2,Form(" B = %f\n", fB));
228 AliDebug(2,Form(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold));
229 AliDebug(2,Form(" Sampling = %f\n", fSampling));
230 AliDebug(2,Form("---------------------------------------------------\n"));
231
814ad4bf 232 return ;
233 }
6569f329 234
235
814ad4bf 236 if(strstr(option,"tim"))
237 gBenchmark->Start("EMCALSDigitizer");
97f3344e 238
33c3c91a 239 AliRunLoader *rl = AliRunLoader::Instance();
5dee926e 240 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
13e1db42 241
88cb7938 242 //switch off reloading of this task while getting event
243 if (!fInit) { // to prevent overwrite existing file
13e1db42 244 AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
88cb7938 245 return ;
fdebddeb 246 }
247
5082063d 248 if (fLastEvent == -1)
5dee926e 249 fLastEvent = rl->GetNumberOfEvents() - 1 ;
1963b290 250 else {
5dee926e 251 fLastEvent = TMath::Min(fLastEvent, rl->GetNumberOfEvents()-1);
1963b290 252 }
5082063d 253 Int_t nEvents = fLastEvent - fFirstEvent + 1;
254
1963b290 255 Int_t ievent;
256 Float_t energy=0.; // de * fSampling - 23-nov-04
5dee926e 257 rl->LoadKinematics();
258 rl->LoadHits("EMCAL");
259
cf24bdc4 260 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
5dee926e 261 rl->GetEvent(ievent);
262 TTree * treeS = emcalLoader->TreeS();
263 if ( !treeS ) {
264 emcalLoader->MakeSDigitsContainer();
265 treeS = emcalLoader->TreeS();
266 }
267 TClonesArray * hits = emcalLoader->Hits() ;
268 TClonesArray * sdigits = emcalLoader->SDigits() ;
814ad4bf 269 sdigits->Clear();
62741b41 270
eb73853b 271 Int_t nSdigits = 0 ;
272 Int_t i;
6569f329 273
274 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
eb73853b 275 for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
276 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
277 AliEMCALDigit * curSDigit = 0 ;
278 AliEMCALDigit * sdigit = 0 ;
279 Bool_t newsdigit = kTRUE;
280
281 // hit->GetId() - Absolute Id number EMCAL segment
282 if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId())
283 energy = hit->GetEnergy() * fSampling; // 23-nov-04
284 if(energy > fECPrimThreshold )
285 // Assign primary number only if deposited energy is significant
286 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
287 hit->GetIparent(), hit->GetId(),
af5bdd85 288 Digitize(energy), hit->GetTime(),
289 -1, energy ) ;
62741b41 290 else
291 curSDigit = new AliEMCALDigit( -1 ,
292 -1 ,
293 hit->GetId(),
af5bdd85 294 Digitize(energy), hit->GetTime(),
295 -1, energy ) ;
eb73853b 296 } else {
297 Warning("Exec"," abs id %i is bad \n", hit->GetId());
298 newsdigit = kFALSE;
299 curSDigit = 0;
300 }
301
302 if(curSDigit != 0){
303 for(Int_t check= 0; check < nSdigits ; check++) {
304 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
305
306 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?
307 *sdigit = *sdigit + *curSDigit;
308 newsdigit = kFALSE;
afa0c708 309 }
afa0c708 310 }
eb73853b 311 }
312 if (newsdigit) {
313 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
314 nSdigits++ ;
315 }
316 delete curSDigit ;
317 } // loop over all hits (hit = deposited energy/entering particle)
62741b41 318 sdigits->Sort() ;
319
320 nSdigits = sdigits->GetEntriesFast() ;
321 fSDigitsInRun += nSdigits ;
ab6a174f 322
62741b41 323 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
324 AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
325 sdigit->SetIndexInList(i) ;
326 }
327
88cb7938 328 // Now write SDigits
eb73853b 329
62741b41 330 Int_t bufferSize = 32000 ;
5dee926e 331 TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
332 if (sdigitsBranch)
333 sdigitsBranch->SetAddress(&sdigits);
334 else
335 treeS->Branch("EMCAL",&sdigits,bufferSize);
eb73853b 336
5dee926e 337 treeS->Fill();
338
339 emcalLoader->WriteSDigits("OVERWRITE");
62741b41 340
341 //NEXT - SDigitizer
5dee926e 342 emcalLoader->WriteSDigitizer("OVERWRITE"); // why in event cycle ?
62741b41 343
344 if(strstr(option,"deb"))
345 PrintSDigits(option) ;
1963b290 346 }
eb73853b 347
88cb7938 348 Unload();
349
5dee926e 350 emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
88cb7938 351
61e0abb5 352 if(strstr(option,"tim")){
9859bfc0 353 gBenchmark->Stop("EMCALSDigitizer");
1963b290 354 printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n",
eb73853b 355 gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ;
62741b41 356 }
61e0abb5 357}
05a92d59 358
bda7cf56 359//__________________________________________________________________
bda7cf56 360Int_t AliEMCALSDigitizer::Digitize(Float_t energy)const {
361 // Digitize the energy
ab6a174f 362 //
363 //JLK 26-Jun-2008 EXPLANATION LONG OVERDUE:
364 //
365 //We have to digitize the SDigit energy so that it can be stored in
366 //AliEMCALDigit, which has only an ADC amplitude member and
367 //(rightly) no energy member. This method converts the energy to an
368 //integer which can be re-converted back to an energy with the
369 //Calibrate(energy) method when it is time to create Digits from SDigits
370 //
371 Double_t aSignal = fA + energy*fB;
372 if (TMath::Abs(aSignal)>2147483647.0) {
373 //PH 2147483647 is the max. integer
374 //PH This apparently is a problem which needs investigation
375 AliWarning(Form("Too big or too small energy %f",aSignal));
376 aSignal = TMath::Sign((Double_t)2147483647,aSignal);
bda7cf56 377 }
ab6a174f 378 return (Int_t ) aSignal;
379}
bda7cf56 380
381//__________________________________________________________________
ab6a174f 382Float_t AliEMCALSDigitizer::Calibrate(Int_t amp)const {
383 //
384 // Convert the amplitude back to energy in GeV
385 //
386 //JLK 26-Jun-2008 EXPLANATION LONG OVERDUE:
387 //
388 //We have to digitize the SDigit energy with the method Digitize()
389 //so that it can be stored in AliEMCALDigit, which has only an ADC
390 //amplitude member and (rightly) no energy member. This method is
391 //just the reverse of Digitize(): it converts the stored amplitude
392 //back to an energy value in GeV so that the SDigit energies can be
393 //summed before adding noise and creating digits out of them
394 //
395 return (Float_t)(amp - fA)/fB;
396
397}
398
106fc2fa 399
ab6a174f 400//__________________________________________________________________
1963b290 401void AliEMCALSDigitizer::Print1(Option_t * option)
402{
403 Print();
404 PrintSDigits(option);
405}
406
61e0abb5 407//__________________________________________________________________
5dee926e 408void AliEMCALSDigitizer::Print(Option_t *option) const
05a92d59 409{
410 // Prints parameters of SDigitizer
e52475ed 411 printf("Print: \n------------------- %s ------------- option %s\n", GetName() , option) ;
1963b290 412 printf(" fInit %i\n", int(fInit));
413 printf(" fFirstEvent %i\n", fFirstEvent);
414 printf(" fLastEvent %i\n", fLastEvent);
88cb7938 415 printf(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()) ;
1963b290 416 printf(" with digitization parameters A = %f\n", fA) ;
417 printf(" B = %f\n", fB) ;
418 printf(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold) ;
419 printf(" Sampling = %f\n", fSampling);
88cb7938 420 printf("---------------------------------------------------\n") ;
61e0abb5 421}
05a92d59 422
61e0abb5 423//__________________________________________________________________
05a92d59 424Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
425{
426 // Equal operator.
427 // SDititizers are equal if their pedestal, slope and threshold are equal
89e103bd 428 if( (fA==sd.fA)&&(fB==sd.fB)&&
fdebddeb 429 (fECPrimThreshold==sd.fECPrimThreshold))
61e0abb5 430 return kTRUE ;
431 else
432 return kFALSE ;
433}
173558f2 434
61e0abb5 435//__________________________________________________________________
d64c959b 436void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
437{
61e0abb5 438 //Prints list of digits produced at the current pass of AliEMCALDigitizer
7ea6391b 439
33c3c91a 440 AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
5dee926e 441 const TClonesArray * sdigits = rl->SDigits() ;
61e0abb5 442
fdebddeb 443 printf("\n") ;
5dee926e 444 printf("event %i", rl->GetRunLoader()->GetEventNumber());
1963b290 445 printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast());
814ad4bf 446 if(strstr(option,"all")||strstr(option,"EMC")){
9859bfc0 447
61e0abb5 448 //loop over digits
449 AliEMCALDigit * digit;
fdebddeb 450 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1963b290 451 Int_t index, isum=0;
11f9c5ff 452 char * tempo = new char[8192];
814ad4bf 453 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
05a92d59 454 digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
11f9c5ff 455 sprintf(tempo, "\n%6d %8d %6.5e %4d %2d :",
456 digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
d804d556 457 printf("%s",tempo);
1963b290 458 isum += digit->GetAmp();
61e0abb5 459
460 Int_t iprimary;
9859bfc0 461 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
d804d556 462 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
463 printf("%s",tempo);
9859bfc0 464 }
61e0abb5 465 }
11f9c5ff 466 delete tempo ;
ab6a174f 467 printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, Calibrate(isum));
1963b290 468 } else printf("\n");
61e0abb5 469}
05a92d59 470
05a92d59 471//____________________________________________________________________________
88cb7938 472void AliEMCALSDigitizer::Unload() const
05a92d59 473{
d64c959b 474 // Unload Hits and SDigits from the folder
33c3c91a 475 AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
5dee926e 476 rl->UnloadHits() ;
477 rl->UnloadSDigits() ;
05a92d59 478}
1963b290 479
fe93d365 480//____________________________________________________________________________
1963b290 481void AliEMCALSDigitizer::Browse(TBrowser* b)
482{
1963b290 483 TTask::Browse(b);
484}