remove redundant QA classes
[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 ---
bda7cf56 52#include <TBenchmark.h>
53#include <TH1.h>
54#include <TBrowser.h>
55#include <Riostream.h>
56#include <TMath.h>
a1e17193 57#include <TROOT.h>
05a92d59 58
61e0abb5 59// --- Standard library ---
d64c959b 60#include "stdlib.h"
61e0abb5 61
62// --- AliRoot header files ---
13e1db42 63#include "AliLog.h"
5dee926e 64#include "AliRunLoader.h"
65#include "AliStack.h"
61e0abb5 66#include "AliEMCALDigit.h"
5dee926e 67#include "AliEMCALLoader.h"
05a92d59 68#include "AliEMCALHit.h"
69#include "AliEMCALSDigitizer.h"
88cb7938 70#include "AliEMCALGeometry.h"
7ea6391b 71//JLK
72//#include "AliEMCALHistoUtilities.h"
61e0abb5 73
74ClassImp(AliEMCALSDigitizer)
61e0abb5 75
76//____________________________________________________________________________
18a21c7c 77AliEMCALSDigitizer::AliEMCALSDigitizer()
78 : TTask("",""),
79 fA(0.),fB(0.),fECPrimThreshold(0.),
80 fDefaultInit(kTRUE),
81 fEventFolderName(0),
82 fInit(0),
83 fSDigitsInRun(0),
84 fFirstEvent(0),
85 fLastEvent(0),
7ea6391b 86 fSampling(0.)
87 //JLK
88 //fControlHists(0),
89 //fHists(0)
61e0abb5 90{
91 // ctor
0a4cb131 92 InitParameters();
61e0abb5 93}
94
95//____________________________________________________________________________
e191bb57 96AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName,
18a21c7c 97 const char * eventFolderName)
98 : TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
99 fA(0.),fB(0.),fECPrimThreshold(0.),
100 fDefaultInit(kFALSE),
101 fEventFolderName(eventFolderName),
102 fInit(0),
103 fSDigitsInRun(0),
104 fFirstEvent(0),
105 fLastEvent(0),
7ea6391b 106 fSampling(0.)
107 //JLK
108 //fControlHists(1),
109 //fHists(0)
61e0abb5 110{
111 // ctor
814ad4bf 112 Init();
62741b41 113 InitParameters() ;
7ea6391b 114
115 //JLK
116 //if(fControlHists) BookControlHists(1);
61e0abb5 117}
118
88cb7938 119
61e0abb5 120//____________________________________________________________________________
18a21c7c 121AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd)
122 : TTask(sd.GetName(),sd.GetTitle()),
123 fA(sd.fA),
124 fB(sd.fB),
125 fECPrimThreshold(sd.fECPrimThreshold),
126 fDefaultInit(sd.fDefaultInit),
127 fEventFolderName(sd.fEventFolderName),
128 fInit(sd.fInit),
129 fSDigitsInRun(sd.fSDigitsInRun),
130 fFirstEvent(sd.fFirstEvent),
131 fLastEvent(sd.fLastEvent),
7ea6391b 132 fSampling(sd.fSampling)
133 //JLK
134 //fControlHists(sd.fControlHists),
135 //fHists(sd.fHists)
18a21c7c 136{
88cb7938 137 //cpy ctor
61e0abb5 138}
36657114 139
88cb7938 140
61e0abb5 141//____________________________________________________________________________
797e311f 142AliEMCALSDigitizer::~AliEMCALSDigitizer() {
14ce0a6e 143 //dtor
5dee926e 144 AliLoader *emcalLoader=0;
145 if ((emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")))
146 emcalLoader->CleanSDigitizer();
797e311f 147}
148
149//____________________________________________________________________________
61e0abb5 150void AliEMCALSDigitizer::Init(){
814ad4bf 151 // Initialization: open root-file, allocate arrays for hits and sdigits,
36657114 152 // attach task SDigitizer to the list of EMCAL tasks
814ad4bf 153 //
154 // Initialization can not be done in the default constructor
155 //============================================================= YS
156 // The initialisation is now done by the getter
61e0abb5 157
dc7da436 158 fInit = kTRUE;
afa0c708 159
5dee926e 160 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
161
162 if ( emcalLoader == 0 ) {
163 Fatal("Init", "Could not obtain the AliEMCALLoader");
814ad4bf 164 return ;
165 }
166
5dee926e 167 emcalLoader->PostSDigitizer(this);
168 emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
88cb7938 169
839828a6 170}
171
172//____________________________________________________________________________
afa0c708 173void AliEMCALSDigitizer::InitParameters()
174{
14ce0a6e 175 //initialize parameters for sdigitization
176
5dee926e 177 const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
62741b41 178 if (geom->GetSampling() == 0.) {
5082063d 179 Fatal("InitParameters", "Sampling factor not set !") ;
62741b41 180 }
1963b290 181 // Initializes parameters
182 fA = 0;
63b36084 183 fB = 1.e+6; // Changed 24 Apr 2007. Dynamic range now 2 TeV
1963b290 184 fSampling = geom->GetSampling();
185
ff98b189 186 // threshold for deposit energy of hit
187 fECPrimThreshold = 0.05;// GeV // 22-may-07 was 0// 24-nov-04 - was 1.e-6;
188
b42d4197 189 AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
190 AliDebug(2,Form(" fInit %i\n", int(fInit)));
191 AliDebug(2,Form(" fFirstEvent %i\n", fFirstEvent));
192 AliDebug(2,Form(" fLastEvent %i\n", fLastEvent));
193 AliDebug(2,Form(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()));
194 AliDebug(2,Form(" with digitization parameters A = %f\n", fA));
195 AliDebug(2,Form(" B = %f\n", fB));
196 AliDebug(2,Form(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold));
197 AliDebug(2,Form(" Sampling = %f\n", fSampling));
198 AliDebug(2,Form("---------------------------------------------------\n"));
199
200 // Print("");
201
839828a6 202}
203
61e0abb5 204//____________________________________________________________________________
afa0c708 205void AliEMCALSDigitizer::Exec(Option_t *option)
206{
fdebddeb 207 // Collects all hit of the same tower into digits
1963b290 208 TString o(option); o.ToUpper();
814ad4bf 209 if (strstr(option, "print") ) {
b42d4197 210
211 AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
212 AliDebug(2,Form(" fInit %i\n", int(fInit)));
213 AliDebug(2,Form(" fFirstEvent %i\n", fFirstEvent));
214 AliDebug(2,Form(" fLastEvent %i\n", fLastEvent));
215 AliDebug(2,Form(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()));
216 AliDebug(2,Form(" with digitization parameters A = %f\n", fA));
217 AliDebug(2,Form(" B = %f\n", fB));
218 AliDebug(2,Form(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold));
219 AliDebug(2,Form(" Sampling = %f\n", fSampling));
220 AliDebug(2,Form("---------------------------------------------------\n"));
221
222 // Print();
814ad4bf 223 return ;
224 }
61e0abb5 225
814ad4bf 226 if(strstr(option,"tim"))
227 gBenchmark->Start("EMCALSDigitizer");
97f3344e 228
5dee926e 229 AliRunLoader *rl = AliRunLoader::GetRunLoader();
230 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
13e1db42 231
88cb7938 232 //switch off reloading of this task while getting event
233 if (!fInit) { // to prevent overwrite existing file
13e1db42 234 AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
88cb7938 235 return ;
fdebddeb 236 }
237
5082063d 238 if (fLastEvent == -1)
5dee926e 239 fLastEvent = rl->GetNumberOfEvents() - 1 ;
1963b290 240 else {
5dee926e 241 fLastEvent = TMath::Min(fLastEvent, rl->GetNumberOfEvents()-1);
1963b290 242 }
5082063d 243 Int_t nEvents = fLastEvent - fFirstEvent + 1;
244
1963b290 245 Int_t ievent;
246 Float_t energy=0.; // de * fSampling - 23-nov-04
5dee926e 247 rl->LoadKinematics();
248 rl->LoadHits("EMCAL");
249
cf24bdc4 250 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
5dee926e 251 rl->GetEvent(ievent);
252 TTree * treeS = emcalLoader->TreeS();
253 if ( !treeS ) {
254 emcalLoader->MakeSDigitsContainer();
255 treeS = emcalLoader->TreeS();
256 }
257 TClonesArray * hits = emcalLoader->Hits() ;
258 TClonesArray * sdigits = emcalLoader->SDigits() ;
814ad4bf 259 sdigits->Clear();
62741b41 260
eb73853b 261 Int_t nSdigits = 0 ;
262 Int_t i;
263 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
264 for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
265 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
266 AliEMCALDigit * curSDigit = 0 ;
267 AliEMCALDigit * sdigit = 0 ;
268 Bool_t newsdigit = kTRUE;
269
270 // hit->GetId() - Absolute Id number EMCAL segment
271 if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId())
272 energy = hit->GetEnergy() * fSampling; // 23-nov-04
273 if(energy > fECPrimThreshold )
274 // Assign primary number only if deposited energy is significant
275 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
276 hit->GetIparent(), hit->GetId(),
af5bdd85 277 Digitize(energy), hit->GetTime(),
278 -1, energy ) ;
62741b41 279 else
280 curSDigit = new AliEMCALDigit( -1 ,
281 -1 ,
282 hit->GetId(),
af5bdd85 283 Digitize(energy), hit->GetTime(),
284 -1, energy ) ;
eb73853b 285 } else {
286 Warning("Exec"," abs id %i is bad \n", hit->GetId());
287 newsdigit = kFALSE;
288 curSDigit = 0;
289 }
290
291 if(curSDigit != 0){
292 for(Int_t check= 0; check < nSdigits ; check++) {
293 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
294
295 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?
296 *sdigit = *sdigit + *curSDigit;
297 newsdigit = kFALSE;
afa0c708 298 }
afa0c708 299 }
eb73853b 300 }
301 if (newsdigit) {
302 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
303 nSdigits++ ;
304 }
305 delete curSDigit ;
306 } // loop over all hits (hit = deposited energy/entering particle)
62741b41 307 sdigits->Sort() ;
308
309 nSdigits = sdigits->GetEntriesFast() ;
310 fSDigitsInRun += nSdigits ;
eb73853b 311
7ea6391b 312 //JLK
313 //Double_t e=0.,esum=0.;
314 //AliEMCALHistoUtilities::FillH1(fHists, 0, double(sdigits->GetEntriesFast()));
62741b41 315 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
316 AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
317 sdigit->SetIndexInList(i) ;
eb73853b 318
7ea6391b 319 //JLK
320 //AliEMCALHistoUtilities::FillH1(fHists, 2, double(sdigit->GetAmp()));
321 //e = double(Calibrate(sdigit->GetAmp()));
322 //esum += e;
323 //AliEMCALHistoUtilities::FillH1(fHists, 3, e);
324 //AliEMCALHistoUtilities::FillH1(fHists, 4, double(sdigit->GetId()));
62741b41 325 }
7ea6391b 326 //JLK
327 //if(esum>0.) AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
62741b41 328
88cb7938 329 // Now write SDigits
eb73853b 330
62741b41 331 Int_t bufferSize = 32000 ;
5dee926e 332 TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
333 if (sdigitsBranch)
334 sdigitsBranch->SetAddress(&sdigits);
335 else
336 treeS->Branch("EMCAL",&sdigits,bufferSize);
eb73853b 337
5dee926e 338 treeS->Fill();
339
340 emcalLoader->WriteSDigits("OVERWRITE");
62741b41 341
342 //NEXT - SDigitizer
5dee926e 343 emcalLoader->WriteSDigitizer("OVERWRITE"); // why in event cycle ?
62741b41 344
345 if(strstr(option,"deb"))
346 PrintSDigits(option) ;
1963b290 347 }
eb73853b 348
88cb7938 349 Unload();
350
5dee926e 351 emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
88cb7938 352
61e0abb5 353 if(strstr(option,"tim")){
9859bfc0 354 gBenchmark->Stop("EMCALSDigitizer");
1963b290 355 printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n",
eb73853b 356 gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ;
62741b41 357 }
61e0abb5 358}
05a92d59 359
bda7cf56 360//__________________________________________________________________
bda7cf56 361Int_t AliEMCALSDigitizer::Digitize(Float_t energy)const {
362 // Digitize the energy
363 Double_t aSignal = fA + energy*fB;
364 if (TMath::Abs(aSignal)>2147483647.0) {
365 //PH 2147483647 is the max. integer
366 //PH This apparently is a problem which needs investigation
367 AliWarning(Form("Too big or too small energy %f",aSignal));
368 aSignal = TMath::Sign((Double_t)2147483647,aSignal);
369 }
370 return (Int_t ) aSignal;
371 }
372
373
374//__________________________________________________________________
106fc2fa 375
1963b290 376void AliEMCALSDigitizer::Print1(Option_t * option)
377{
378 Print();
379 PrintSDigits(option);
380}
381
61e0abb5 382//__________________________________________________________________
5dee926e 383void AliEMCALSDigitizer::Print(Option_t *option) const
05a92d59 384{
385 // Prints parameters of SDigitizer
e52475ed 386 printf("Print: \n------------------- %s ------------- option %s\n", GetName() , option) ;
1963b290 387 printf(" fInit %i\n", int(fInit));
388 printf(" fFirstEvent %i\n", fFirstEvent);
389 printf(" fLastEvent %i\n", fLastEvent);
88cb7938 390 printf(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()) ;
1963b290 391 printf(" with digitization parameters A = %f\n", fA) ;
392 printf(" B = %f\n", fB) ;
393 printf(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold) ;
394 printf(" Sampling = %f\n", fSampling);
88cb7938 395 printf("---------------------------------------------------\n") ;
61e0abb5 396}
05a92d59 397
61e0abb5 398//__________________________________________________________________
05a92d59 399Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
400{
401 // Equal operator.
402 // SDititizers are equal if their pedestal, slope and threshold are equal
89e103bd 403 if( (fA==sd.fA)&&(fB==sd.fB)&&
fdebddeb 404 (fECPrimThreshold==sd.fECPrimThreshold))
61e0abb5 405 return kTRUE ;
406 else
407 return kFALSE ;
408}
173558f2 409
61e0abb5 410//__________________________________________________________________
d64c959b 411void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
412{
61e0abb5 413 //Prints list of digits produced at the current pass of AliEMCALDigitizer
7ea6391b 414
5dee926e 415 AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
416 const TClonesArray * sdigits = rl->SDigits() ;
61e0abb5 417
fdebddeb 418 printf("\n") ;
5dee926e 419 printf("event %i", rl->GetRunLoader()->GetEventNumber());
1963b290 420 printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast());
814ad4bf 421 if(strstr(option,"all")||strstr(option,"EMC")){
9859bfc0 422
61e0abb5 423 //loop over digits
424 AliEMCALDigit * digit;
fdebddeb 425 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1963b290 426 Int_t index, isum=0;
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()) ;
1963b290 432 printf(tempo);
433 isum += digit->GetAmp();
61e0abb5 434
435 Int_t iprimary;
9859bfc0 436 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
11f9c5ff 437 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
fdebddeb 438 printf(tempo);
9859bfc0 439 }
61e0abb5 440 }
11f9c5ff 441 delete tempo ;
1963b290 442 printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6);
443 } else printf("\n");
61e0abb5 444}
05a92d59 445
05a92d59 446//____________________________________________________________________________
88cb7938 447void AliEMCALSDigitizer::Unload() const
05a92d59 448{
d64c959b 449 // Unload Hits and SDigits from the folder
5dee926e 450 AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
451 rl->UnloadHits() ;
452 rl->UnloadSDigits() ;
05a92d59 453}
1963b290 454
fe93d365 455//____________________________________________________________________________
1963b290 456void AliEMCALSDigitizer::Browse(TBrowser* b)
457{
7ea6391b 458 //JLK
459 //if(fHists) b->Add(fHists);
1963b290 460 TTask::Browse(b);
461}
462
7ea6391b 463/*
fe93d365 464//____________________________________________________________________________
eb0b1051 465TList *AliEMCALSDigitizer::BookControlHists(int var)
14ce0a6e 466{
467 //book histograms for monitoring sdigitization
468 // 22-nov-04
1963b290 469 gROOT->cd();
5dee926e 470 const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance() ;
1963b290 471 if(var>=1){
4bba84bd 472 AliDebug(1, " BookControlHists() in action ");
1963b290 473 new TH1F("HSDigiN", "#EMCAL sdigits ", 1001, -0.5, 1000.5);
474 new TH1F("HSDigiSumEnergy","Sum.EMCAL energy", 1000, 0.0, 100.);
475 new TH1F("HSDigiAmp", "EMCAL sdigits amplitude", 1000, 0., 2.e+9);
476 new TH1F("HSDigiEnergy","EMCAL cell energy", 1000, 0.0, 100.);
477 new TH1F("HSDigiAbsId","EMCAL absID for sdigits",
478 geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
479 }
e52475ed 480
315d1c64 481 fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalSDigiControlHists", kFALSE);
4bba84bd 482 // fHists = 0; ??
e52475ed 483
1963b290 484 return fHists;
485}
486
fe93d365 487//____________________________________________________________________________
1963b290 488void AliEMCALSDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
489{
315d1c64 490 AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt);
1963b290 491}
7ea6391b 492*/