Change thresholds ans noise parameters
[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
61e0abb5 46//////////////////////////////////////////////////////////////////////////////
47
48
49// --- ROOT system ---
50#include "TFile.h"
51#include "TTask.h"
52#include "TTree.h"
53#include "TSystem.h"
54#include "TROOT.h"
55#include "TFolder.h"
56#include "TBenchmark.h"
57// --- Standard library ---
58#include <iomanip.h>
59
60// --- AliRoot header files ---
61#include "AliRun.h"
62#include "AliEMCALDigit.h"
63#include "AliEMCALHit.h"
64#include "AliEMCALSDigitizer.h"
65#include "AliEMCALGeometry.h"
66#include "AliEMCALv1.h"
814ad4bf 67#include "AliEMCALGetter.h"
61e0abb5 68
69ClassImp(AliEMCALSDigitizer)
70
71
72//____________________________________________________________________________
73 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
74{
75 // ctor
76 fA = 0;
e1f60236 77 fB = 0. ;
78 fPrimThreshold = 0. ;
61e0abb5 79 fNevents = 0 ;
80 fSDigits = 0 ;
81 fHits = 0 ;
e1f60236 82 fLayerRatio = 0. ;
61e0abb5 83 fIsInitialized = kFALSE ;
84
85}
86
87//____________________________________________________________________________
814ad4bf 88AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask(sDigitsTitle, headerFile)
61e0abb5 89{
90 // ctor
91 fA = 0;
92 fB = 10000000.;
ffa6d63b 93 fPrimThreshold = 0.01 ;
814ad4bf 94 fNevents = 0 ;
e1f60236 95 fLayerRatio = 5./6. ;
814ad4bf 96 Init();
61e0abb5 97}
98
99//____________________________________________________________________________
100AliEMCALSDigitizer::~AliEMCALSDigitizer()
101{
102 // dtor
103 if(fSDigits)
104 delete fSDigits ;
105 if(fHits)
106 delete fHits ;
107}
108//____________________________________________________________________________
109void AliEMCALSDigitizer::Init(){
61e0abb5 110
814ad4bf 111 // Initialization: open root-file, allocate arrays for hits and sdigits,
112 // attach task SDigitizer to the list of PHOS tasks
113 //
114 // Initialization can not be done in the default constructor
115 //============================================================= YS
116 // The initialisation is now done by the getter
61e0abb5 117
814ad4bf 118if( strcmp(GetTitle(), "") == 0 )
119 SetTitle("galice.root") ;
120
121 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), "update") ;
122 if ( gime == 0 ) {
123 cerr << "ERROR: AliEMCALSDigitizer::Init -> Could not obtain the Getter object !" << endl ;
124 return ;
125 }
126
127 gime->PostSDigits( GetName(), GetTitle() ) ;
61e0abb5 128
814ad4bf 129 TString sdname(GetName() );
130 sdname.Append(":") ;
131 sdname.Append(GetTitle() ) ;
132 SetName(sdname) ;
133 gime->PostSDigitizer(this) ;
134
135 // create a folder on the white board //YSAlice/WhiteBoard/SDigits/EMCAL/headerFile/sdigitsTitle
136
137 }
61e0abb5 138//____________________________________________________________________________
139void AliEMCALSDigitizer::Exec(Option_t *option) {
61e0abb5 140
814ad4bf 141// Collects all hits in the same active volume into digit
142
143 if( strcmp(GetName(), "") == 0 )
144 Init() ;
61e0abb5 145
814ad4bf 146 if (strstr(option, "print") ) {
147 Print("") ;
148 return ;
149 }
61e0abb5 150
814ad4bf 151 if(strstr(option,"tim"))
152 gBenchmark->Start("EMCALSDigitizer");
153 //Check, if this branch already exits
154 gAlice->GetEvent(0) ;
155 if(gAlice->TreeS() ) {
156 TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
157 TIter next(lob) ;
158 TBranch * branch = 0 ;
159 Bool_t emcalfound = kFALSE, sdigitizerfound = kFALSE ;
61e0abb5 160
814ad4bf 161 while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
162 if ( (strcmp(branch->GetName(), "EMCAL")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
163 emcalfound = kTRUE ;
164
165 else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
166 sdigitizerfound = kTRUE ;
61e0abb5 167 }
168
814ad4bf 169 if ( emcalfound || sdigitizerfound ) {
170 cerr << "WARNING: AliEMCALSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName()
171 << " already exits" << endl ;
172 return ;
173 }
174 }
175 TString sdname(GetName()) ;
176 sdname.Remove(sdname.Index(GetTitle())-1) ;
177 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
814ad4bf 178 Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
179
180 Int_t ievent ;
181 for(ievent = 0; ievent < nevents; ievent++){
182 gime->Event(ievent,"H") ;
183 const TClonesArray * fHits = gime->Hits() ;
184 TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
185 sdigits->Clear();
61e0abb5 186 Int_t nSdigits = 0 ;
814ad4bf 187
188 //Collects all hits in the same active volume into digit
189
190
61e0abb5 191
192
193 //Now made SDigits from hits, for EMCAL it is the same, so just copy
814ad4bf 194 // Int_t itrack ;
195 // for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
196 //gime->Track(itrack);
61e0abb5 197 //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
814ad4bf 198
199
200
61e0abb5 201
202 Int_t i;
e1f60236 203 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { // loop over all hits (hit = deposited energy/layer/entering particle)
204 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(fHits->At(i)) ;
205 AliEMCALDigit * curSDigit = 0 ;
206 AliEMCALDigit * sdigit = 0 ;
ffa6d63b 207 Bool_t newsdigit = kTRUE;
e1f60236 208
209 // Assign primary number only if deposited energy is significant
210
211 Float_t preShowerFactor ;
212 if (hit->IsInPreShower() )
213 preShowerFactor = fLayerRatio ;
214 else
215 preShowerFactor = 1 ;
216
229f77c4 217 if( hit->GetEnergy() > fPrimThreshold)//{
218 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
e1f60236 219 hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
220 Digitize( hit->GetEnergy() * preShowerFactor ) ,
814ad4bf 221 hit->GetTime()) ;
229f77c4 222 //cout << "Sdigitizer > Threshold Hit Primary = " << hit->GetPrimary() << "Digit Primary = "<< curSDigit->GetPrimary(1) << endl ;}
223 else //{
814ad4bf 224 curSDigit = new AliEMCALDigit( -1 ,
225 -1 ,
e1f60236 226 Layer2TowerID(hit->GetId(),hit->IsInPreShower()),
227 Digitize( hit->GetEnergy() * preShowerFactor ) ,
228 hit->GetTime() ) ;
229f77c4 229 // cout << "SDigitizer < threshold Hit Primary = " << hit->GetPrimary() << "Digit Primary = "<< curSDigit->GetPrimary(1) << endl ;}
e1f60236 230 Int_t check = 0 ;
231 for(check= 0; check < nSdigits ; check++) {
814ad4bf 232 sdigit = (AliEMCALDigit *)sdigits->At(check);
e1f60236 233 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same tower or the same preshower ?
234 *sdigit = *sdigit + *curSDigit;
235
236 newsdigit = kFALSE;
237 }
814ad4bf 238 }
e1f60236 239 if (newsdigit) {
240 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
ffa6d63b 241 nSdigits++ ;
e1f60236 242 }
243 } // loop over all hits (hit = deposited energy/layer/entering particle)
814ad4bf 244 // } // loop over tracks
e1f60236 245
814ad4bf 246 sdigits->Sort() ;
61e0abb5 247
814ad4bf 248 nSdigits = sdigits->GetEntriesFast() ;
249 sdigits->Expand(nSdigits) ;
250
251 // Int_t i ;
e1f60236 252 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
253
229f77c4 254 Int_t lastPreShowerIndex = nSdigits - 1 ;
255 if (!(dynamic_cast<AliEMCALDigit *>(sdigits->At(lastPreShowerIndex))->IsInPreShower()))
256 lastPreShowerIndex = -2;
e1f60236 257 Int_t firstPreShowerIndex = -1 ;
258 Int_t index ;
259 AliEMCALDigit * sdigit = 0 ;
260 for ( index = 0; index < nSdigits ; index++) {
261 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ) ;
262 if (sdigit->IsInPreShower() ){
263 firstPreShowerIndex = index ;
264 break ;
265 }
266 }
267
268 AliEMCALDigit * preshower ;
269 AliEMCALDigit * tower ;
270 Int_t lastIndex = lastPreShowerIndex +1 ;
271
272
273 for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) {
274 preshower = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) );
275 Bool_t towerFound = kFALSE ;
276 Int_t jndex ;
277 for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
278 tower = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) );
279 if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {
280 *tower = *tower + *preshower ; // and add preshower to tower
281 towerFound = kTRUE ;
282 }
283 }
284 if ( !towerFound ) {
285
286 new((*sdigits)[lastIndex]) AliEMCALDigit(*preshower);
287 AliEMCALDigit * temp = dynamic_cast<AliEMCALDigit *>(sdigits->At(lastIndex)) ;
288 temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ;
289 lastIndex++ ;
290 }
291 }
229f77c4 292
e1f60236 293 sdigits->Sort() ;
229f77c4 294 Int_t NPrimarymax = -1 ;
295 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
296 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
297 sdigit->SetIndexInList(i) ;
298 }
299
300 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
301 if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
302 NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
303 }
304 if(gAlice->TreeS() == 0)
814ad4bf 305 gAlice->MakeTree("S") ;
229f77c4 306 cout << " AliEMCALSDigitizer:: NPrimaryMax = " << NPrimarymax << endl;
814ad4bf 307
308
309
310 //Make (if necessary) branches
311 char * file =0;
312 if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
313 file = new char[strlen(gAlice->GetBaseFile())+20] ;
314 sprintf(file,"%s/EMCAL.SDigits.root",gAlice->GetBaseFile()) ;
61e0abb5 315 }
61e0abb5 316
814ad4bf 317 TDirectory *cwd = gDirectory;
318
319 //First list of sdigits
320 Int_t bufferSize = 32000 ;
321 TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
322 sdigitsBranch->SetTitle(sdname);
323 cout << " AliEMCALSDigitizer::Exec sdname " << sdname << endl ;
324
325 if (file) {
326 sdigitsBranch->SetFile(file);
327 TIter next( sdigitsBranch->GetListOfBranches());
328 TBranch * subbr;
329 while ((subbr=(TBranch*)next())) {
330 subbr->SetFile(file);
331 }
332 cwd->cd();
333 }
334
335 //second - SDigitizer
336 Int_t splitlevel = 0 ;
337 AliEMCALSDigitizer * sd = this ;
338 TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
339 &sd,bufferSize,splitlevel);
340 sdigitizerBranch->SetTitle(sdname);
341 if (file) {
342 sdigitizerBranch->SetFile(file);
343 TIter next( sdigitizerBranch->GetListOfBranches());
344 TBranch * subbr ;
345 while ((subbr=(TBranch*)next())) {
346 subbr->SetFile(file);
347 }
348 cwd->cd();
349 delete file;
350 }
351
352 sdigitsBranch->Fill() ;
353 sdigitizerBranch->Fill() ;
354 gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
355
356 if(strstr(option,"deb"))
357 PrintSDigits(option) ;
358
61e0abb5 359 }
360
361 if(strstr(option,"tim")){
362 gBenchmark->Stop("EMCALSDigitizer");
363 cout << "AliEMCALSDigitizer:" << endl ;
364 cout << " took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing "
365 << gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
366 cout << endl ;
367 }
368
369
370}
371//__________________________________________________________________
372void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
814ad4bf 373
374 // Setting title to branch SDigits
375
376 TString stitle(title) ;
377
378 // check if branch with title already exists
379 TBranch * sdigitsBranch =
380 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
381 TBranch * sdigitizerBranch =
382 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
383 const char * sdigitsTitle = sdigitsBranch ->GetTitle() ;
384 const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
385 if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
386 cerr << "ERROR: AliEMCALSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
387 return ;
388 }
389
390 cout << "AliEMCALSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
391
392 SetName(title) ;
393
394 // Post to the WhiteBoard
395 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
396 gime->PostSDigits( title, GetTitle()) ;
397
398
61e0abb5 399}
400//__________________________________________________________________
401void AliEMCALSDigitizer::Print(Option_t* option)const{
402 cout << "------------------- "<< GetName() << " -------------" << endl ;
814ad4bf 403 cout << " Writing SDigitis to branch with title " << GetName() << endl ;
61e0abb5 404 cout << " with digitization parameters A = " << fA << endl ;
405 cout << " B = " << fB << endl ;
406 cout << " Threshold for Primary assignment= " << fPrimThreshold << endl ;
407 cout << "---------------------------------------------------"<<endl ;
408
409}
410//__________________________________________________________________
411Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
412 if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
413 return kTRUE ;
414 else
415 return kFALSE ;
416}
417//__________________________________________________________________
418void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
419 //Prints list of digits produced at the current pass of AliEMCALDigitizer
420
814ad4bf 421 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
422 TString sdname(GetName()) ;
423 sdname.Remove(sdname.Index(GetTitle())-1) ;
424 const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
61e0abb5 425
814ad4bf 426 cout << "AliEMCALSDigitiser: event " << gAlice->GetEvNumber() << endl ;
427 cout << " Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
428 cout << endl ;
429 if(strstr(option,"all")||strstr(option,"EMC")){
430
61e0abb5 431 //loop over digits
432 AliEMCALDigit * digit;
e1f60236 433 cout << "SDigit Id " << " Amplitude " << " Time " << " Index " << " Nprim " << " Primaries list " << endl;
61e0abb5 434 Int_t index ;
814ad4bf 435 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
436 digit = (AliEMCALDigit * ) sdigits->At(index) ;
437 cout << setw(6) << digit->GetId() << " " << setw(10) << digit->GetAmp() << " " << digit->GetTime()
438 << setw(6) << digit->GetIndexInList() << " "
439 << setw(5) << digit->GetNprimary() <<" ";
61e0abb5 440
441 Int_t iprimary;
442 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
814ad4bf 443 cout << " " << digit->GetPrimary(iprimary+1) << " ";
61e0abb5 444 cout << endl;
445 }
814ad4bf 446 cout <<endl;
61e0abb5 447 }
448}
814ad4bf 449//________________________________________________________________________
229f77c4 450Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower){
451 // Method to Transform from Hit Id to Digit Id
452 // This function should be one to one
814ad4bf 453 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
454 const AliEMCALGeometry * geom = gime->EMCALGeometry();
229f77c4 455 Int_t ieta = (ihit/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index
456 Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index
457 Int_t it = -10;
458 Int_t ipre = 0;
814ad4bf 459
229f77c4 460 if (preshower)ipre = 1;
461 if (iphi > 0 && ieta >= 0){
462 it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ();
463 return it;
464 }else{
465 cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Eta number = "
466 << ieta << "Phi number = " << iphi << endl ;
467 return it;
468 } // end if iphi>0 && ieta>0
814ad4bf 469}
229f77c4 470//_______________________________________________________________________________________
471void AliEMCALSDigitizer::TestTowerID(void)
472{
473 Int_t j;
474
475 Bool_t preshower = kFALSE;
476 for (j = 0 ; j < 10 ; j++){ // loop over hit id
477 Int_t i;
478 for (i = 0 ; i <= 2 ; i++){ // loop over
479 Int_t k = i*96*144+j*144+1;
480 cout << " Hit Index = " << k << " " << j*10 << " TOWERID = " << Layer2TowerID(k, preshower) << endl ;
481 }
482 }
814ad4bf 483}