Bug in title-name corrected
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.cxx
CommitLineData
61e0abb5 1/**************************************************************************
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
21// from one primary in one active cell
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//
44//*-- Author : Dmitri Peressounko (SUBATECH & KI)
45//////////////////////////////////////////////////////////////////////////////
46
47
48// --- ROOT system ---
49#include "TFile.h"
50#include "TTask.h"
51#include "TTree.h"
52#include "TSystem.h"
53#include "TROOT.h"
54#include "TFolder.h"
55#include "TBenchmark.h"
56// --- Standard library ---
57#include <iomanip.h>
58
59// --- AliRoot header files ---
60#include "AliRun.h"
61#include "AliEMCALDigit.h"
62#include "AliEMCALHit.h"
63#include "AliEMCALSDigitizer.h"
64#include "AliEMCALGeometry.h"
65#include "AliEMCALv1.h"
66
67ClassImp(AliEMCALSDigitizer)
68
69
70//____________________________________________________________________________
71 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
72{
73 // ctor
74 fA = 0;
75 fB = 10000000. ;
76 fPrimThreshold = 0.001 ;
77 fNevents = 0 ;
78 fSDigits = 0 ;
79 fHits = 0 ;
80 fIsInitialized = kFALSE ;
81
82}
83
84//____________________________________________________________________________
85AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliEMCALSDigitizer","")
86{
87 // ctor
88 fA = 0;
89 fB = 10000000.;
90 fPrimThreshold = 0.001 ;
91 fNevents = 0 ;
92 fSDigitsTitle = sDigitsTitle ;
93 fHeadersFile = headerFile ;
94 fSDigits = new TClonesArray("AliEMCALDigit",1000);
95 fHits = new TClonesArray("AliEMCALHit",1000);
96
97 TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
98
99 //File was not opened yet
100 if(file == 0){
101 if(fHeadersFile.Contains("rfio"))
102 file = TFile::Open(fHeadersFile,"update") ;
103 else
104 file = new TFile(fHeadersFile.Data(),"update") ;
105 gAlice = (AliRun *) file->Get("gAlice") ;
106 }
107
108 //add Task to //root/Tasks folder
109 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
110 roottasks->Add(this) ;
111
112 fIsInitialized = kTRUE ;
113}
114
115//____________________________________________________________________________
116AliEMCALSDigitizer::~AliEMCALSDigitizer()
117{
118 // dtor
119 if(fSDigits)
120 delete fSDigits ;
121 if(fHits)
122 delete fHits ;
123}
124//____________________________________________________________________________
125void AliEMCALSDigitizer::Init(){
126 //Initialization can not be done in the default constructor
127
128 if(!fIsInitialized){
129
130 if(fHeadersFile.IsNull())
131 fHeadersFile="galice.root" ;
132
133 TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
134
135 //if file was not opened yet, read gAlice
136 if(file == 0){
137 file = new TFile(fHeadersFile.Data(),"update") ;
138 gAlice = (AliRun *) file->Get("gAlice") ;
139 }
140
141 fHits = new TClonesArray("AliEMCALHit",1000);
142 fSDigits = new TClonesArray("AliEMCALDigit",1000);
143
144 // add Task to //root/Tasks folder
145 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
146 roottasks->Add(this) ;
147
148 fIsInitialized = kTRUE ;
149 }
150}
151//____________________________________________________________________________
152void AliEMCALSDigitizer::Exec(Option_t *option) {
153 //Collects all hits in the same active volume into digit
154
155 if(!fIsInitialized)
156 Init() ;
157
158 if(strstr(option,"tim"))
159 gBenchmark->Start("EMCALSDigitizer");
160
161 fNevents = (Int_t) gAlice->TreeE()->GetEntries() ;
162
163 Int_t ievent ;
164 for(ievent = 0; ievent < fNevents; ievent++){
165 gAlice->GetEvent(ievent) ;
166 gAlice->SetEvent(ievent) ;
167
168 if(gAlice->TreeH()==0){
169 cout << "AliEMCALSDigitizer: There is no Hit Tree" << endl;
170 return ;
171 }
172
173 //set address of the hits
174 TBranch * branch = gAlice->TreeH()->GetBranch("EMCAL");
175 if (branch)
176 branch->SetAddress(&fHits);
177 else{
178 cout << "ERROR in AliEMCALSDigitizer: "<< endl ;
179 cout << " no branch EMCAL in TreeH"<< endl ;
180 cout << " do nothing " << endl ;
181 return ;
182 }
183
184 fSDigits->Clear();
185 Int_t nSdigits = 0 ;
186
187
188 //Now made SDigits from hits, for EMCAL it is the same, so just copy
189 Int_t itrack ;
190 for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
191
192 //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
193 branch->GetEntry(itrack,0);
194 AliEMCAL * EMCAL = (AliEMCAL *) gAlice->GetDetector("EMCAL") ;
195 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance( EMCAL->GetGeometry()->GetName(), EMCAL->GetGeometry()->GetTitle() );
196
197 Int_t i;
198 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
199 AliEMCALHit * hit = (AliEMCALHit*)fHits->At(i) ;
200 AliEMCALDigit curSDigit = AliEMCALDigit(1,1,1,1);
201 AliEMCALDigit *sdigit ;
202
203 // Assign primary number only if contribution is significant
204
205 if( hit->GetEnergy() > fPrimThreshold)
206 curSDigit = AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
207 else
208 curSDigit = AliEMCALDigit( -1 , -1 ,(((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
209
210 for(Int_t check= 0; check < nSdigits; check++) {
211 sdigit = (AliEMCALDigit *)fSDigits->At(check);
212 if( (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()) == ((sdigit->GetId()/geom->GetNPhi())%geom->GetNZ())) && ((hit->GetId()%geom->GetNPhi()) == (sdigit->GetId()%geom->GetNPhi())))
213 {
214 *sdigit = *sdigit + curSDigit ;
215 }
216 else
217 { new((*fSDigits)[nSdigits]) AliEMCALDigit(curSDigit);
218 nSdigits++ ; }
219 }
220
221
222 if( hit->GetEnergy() > fPrimThreshold)
223 curSDigit = AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), ((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
224 else
225 curSDigit = AliEMCALDigit( -1 , -1 ,((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ;
226
227 if((hit->GetId()/geom->GetNPhi()) < (2*geom->GetNZ()))
228 {
229 for(Int_t check= 0; check < nSdigits; check++) {
230 sdigit = (AliEMCALDigit *)fSDigits->At(check);
231 if( (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()) == ((sdigit->GetId()/geom->GetNPhi())%geom->GetNZ())) && ((hit->GetId()%geom->GetNPhi()) == (sdigit->GetId()%geom->GetNPhi())))
232 {
233 *sdigit = *sdigit + curSDigit ;
234 }
235 else
236 { new((*fSDigits)[nSdigits]) AliEMCALDigit(curSDigit);
237 nSdigits++ ; }
238 }
239 }
240 }
241 } // loop over tracks
242
243 fSDigits->Sort() ;
244
245 nSdigits = fSDigits->GetEntriesFast() ;
246 fSDigits->Expand(nSdigits) ;
247
248 Int_t i ;
249 for (i = 0 ; i < nSdigits ; i++) {
250 AliEMCALDigit * digit = (AliEMCALDigit *) fSDigits->At(i) ;
251 digit->SetIndexInList(i) ;
252 }
253
254 if(gAlice->TreeS() == 0)
255 gAlice->MakeTree("S") ;
256
257 //check, if this branch already exits?
258 TBranch * sdigitsBranch = 0;
259 TBranch * sdigitizerBranch = 0;
260
261 TObjArray * branches = gAlice->TreeS()->GetListOfBranches() ;
262 Int_t ibranch;
263 Bool_t emcalNotFound = kTRUE ;
264 Bool_t sdigitizerNotFound = kTRUE ;
265
266 for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
267
268 if(emcalNotFound){
269 sdigitsBranch=(TBranch *) branches->At(ibranch) ;
270 if( (strcmp("EMCAL",sdigitsBranch->GetName())==0 ) &&
271 (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) )
272 emcalNotFound = kFALSE ;
273 }
274 if(sdigitizerNotFound){
275 sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
276 if( (strcmp(sdigitizerBranch->GetName(),"AliEMCALSDigitizer") == 0)&&
277 (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) )
278 sdigitizerNotFound = kFALSE ;
279 }
280 }
281
282 if(!(sdigitizerNotFound && emcalNotFound)){
283 cout << "AliEMCALSdigitizer error:" << endl ;
284 cout << "Can not overwrite existing branches: do not write" << endl ;
285 return ;
286 }
287
288 //Make (if necessary) branches
289 char * file =0;
290 if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
291 file = new char[strlen(gAlice->GetBaseFile())+20] ;
292 sprintf(file,"%s/EMCAL.SDigits.root",gAlice->GetBaseFile()) ;
293 }
294
295 TDirectory *cwd = gDirectory;
296
297 //First list of sdigits
298 Int_t bufferSize = 32000 ;
299 sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&fSDigits,bufferSize);
300 sdigitsBranch->SetTitle(fSDigitsTitle.Data());
301 if (file) {
302 sdigitsBranch->SetFile(file);
303 TIter next( sdigitsBranch->GetListOfBranches());
304 TBranch * subbr;
305 while ((subbr=(TBranch*)next())) {
306 subbr->SetFile(file);
307 }
308 cwd->cd();
309 }
310
311 //second - SDigitizer
312 Int_t splitlevel = 0 ;
313 AliEMCALSDigitizer * sd = this ;
314 sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
315 &sd,bufferSize,splitlevel);
316 sdigitizerBranch->SetTitle(fSDigitsTitle.Data());
317 if (file) {
318 sdigitizerBranch->SetFile(file);
319 TIter next( sdigitizerBranch->GetListOfBranches());
320 TBranch * subbr ;
321 while ((subbr=(TBranch*)next())) {
322 subbr->SetFile(file);
323 }
324 cwd->cd();
325 delete file;
326 }
327
328 sdigitsBranch->Fill() ;
329 sdigitizerBranch->Fill() ;
330 gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
331
332 if(strstr(option,"deb"))
333 PrintSDigits(option) ;
334
335 }
336
337 if(strstr(option,"tim")){
338 gBenchmark->Stop("EMCALSDigitizer");
339 cout << "AliEMCALSDigitizer:" << endl ;
340 cout << " took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing "
341 << gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
342 cout << endl ;
343 }
344
345
346}
347//__________________________________________________________________
348void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
349 //Seting title to branch SDigits
350 if(!fSDigitsTitle.IsNull())
351 cout << "AliEMCALSdigitizer: changing SDigits file from " <<fSDigitsTitle.Data() << " to " << title << endl ;
352 fSDigitsTitle=title ;
353}
354//__________________________________________________________________
355void AliEMCALSDigitizer::Print(Option_t* option)const{
356 cout << "------------------- "<< GetName() << " -------------" << endl ;
357 cout << " Writing SDigitis to branch with title " << fSDigitsTitle.Data() << endl ;
358 cout << " with digitization parameters A = " << fA << endl ;
359 cout << " B = " << fB << endl ;
360 cout << " Threshold for Primary assignment= " << fPrimThreshold << endl ;
361 cout << "---------------------------------------------------"<<endl ;
362
363}
364//__________________________________________________________________
365Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
366 if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
367 return kTRUE ;
368 else
369 return kFALSE ;
370}
371//__________________________________________________________________
372void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
373 //Prints list of digits produced at the current pass of AliEMCALDigitizer
374
375 cout << "AliEMCALSDigitizer: " << endl ;
376 cout << " Number of entries in SDigits list " << fSDigits->GetEntriesFast() << endl ;
377 cout << endl ;
378
379 if(strstr(option,"all")){// print all digits
380
381 //loop over digits
382 AliEMCALDigit * digit;
383 cout << "SDigit Id " << " Amplitude " << " Index " << " Nprim " << " Primaries list " << endl;
384 Int_t index ;
385 for (index = 0 ; index < fSDigits->GetEntries() ; index++) {
386 digit = (AliEMCALDigit * ) fSDigits->At(index) ;
387 cout << setw(8) << digit->GetId() << " " << setw(3) << digit->GetAmp() << " "
388 << setw(6) << digit->GetIndexInList() << " "
389 << setw(5) << digit->GetNprimary() <<" ";
390
391 Int_t iprimary;
392 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
393 cout << setw(5) << digit->GetPrimary(iprimary+1) << " ";
394 cout << endl;
395 }
396
397 }
398}