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