Fixes for cmake
[u/mrichter/AliRoot.git] / PHOS / AliPHOSDigitizer.cxx
CommitLineData
990119d6 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
702ab87e 18/* History of cvs commits:
19 *
a28333c5 20 * $Log: AliPHOSDigitizer.cxx,v $
21 * Revision 1.104 2007/12/18 09:08:18 hristov
22 * Splitting of the QA maker into simulation and reconstruction dependent parts (Yves)
23 *
04236e67 24 * Revision 1.103 2007/11/07 11:25:06 schutz
25 * Comment out the QA checking before starting digitization
26 *
18c3b8c2 27 * Revision 1.102 2007/10/19 18:04:29 schutz
28 * The standalone QA data maker is called from AliSimulation and AliReconstruction outside the event loop; i.e. re-reading the data. The QA data making in the event loop has been commented out.
29 *
c65c502a 30 * Revision 1.101 2007/10/14 21:08:10 schutz
31 * Introduced the checking of QA results from previous step before entering the event loop
32 *
8661738e 33 * Revision 1.100 2007/10/10 09:05:10 schutz
34 * Changing name QualAss to QA
35 *
b8274834 36 * Revision 1.99 2007/09/30 17:08:20 schutz
37 * Introducing the notion of QA data acquisition cycle (needed by online)
38 *
5b188f2f 39 * Revision 1.98 2007/09/26 14:22:17 cvetan
40 * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
41 *
d76c31f4 42 * Revision 1.97 2007/08/07 14:12:03 kharlov
43 * Quality assurance added (Yves Schutz)
44 *
ddd1a39c 45 * Revision 1.96 2007/04/28 10:43:36 policheh
46 * Dead channels simulation: digit energy sets to 0.
47 *
1c6a163f 48 * Revision 1.95 2007/04/10 07:20:52 kharlov
49 * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
50 *
1a197b77 51 * Revision 1.94 2007/02/01 10:34:47 hristov
52 * Removing warnings on Solaris x86
53 *
ad4aeaf4 54 * Revision 1.93 2006/10/17 13:17:01 kharlov
55 * Replace AliInfo by AliDebug
56 *
d6e8d7d3 57 * Revision 1.92 2006/08/28 10:01:56 kharlov
58 * Effective C++ warnings fixed (Timur Pocheptsov)
59 *
3663622c 60 * Revision 1.91 2006/04/29 20:25:30 hristov
61 * Decalibration is implemented (Yu.Kharlov)
62 *
877695e7 63 * Revision 1.90 2006/04/22 10:30:17 hristov
64 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
65 *
27a73a5d 66 * Revision 1.89 2006/04/11 15:22:59 hristov
67 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
68 *
28871337 69 * Revision 1.88 2006/03/13 14:05:43 kharlov
70 * Calibration objects for EMC and CPV
71 *
fc6706cb 72 * Revision 1.87 2005/08/24 15:33:49 kharlov
73 * Calibration data for raw digits
74 *
a8ec0771 75 * Revision 1.86 2005/07/12 20:07:35 hristov
76 * Changes needed to run simulation and reconstrruction in the same AliRoot session
77 *
7c193632 78 * Revision 1.85 2005/05/28 14:19:04 schutz
79 * Compilation warnings fixed by T.P.
80 *
702ab87e 81 */
88cb7938 82
990119d6 83//_________________________________________________________________________
990119d6 84//*-- Author : Dmitri Peressounko (SUBATECH & Kurchatov Institute)
85//////////////////////////////////////////////////////////////////////////////
a4e98857 86// This TTask performs digitization of Summable digits (in the PHOS case it is just
87// the sum of contributions from all primary particles into a given cell).
990119d6 88// In addition it performs mixing of summable digits from different events.
7b7c1533 89// The name of the TTask is also the title of the branch that will contain
90// the created SDigits
91// The title of the TTAsk is the name of the file that contains the hits from
92// which the SDigits are created
bca3b32a 93//
94// For each event two branches are created in TreeD:
95// "PHOS" - list of digits
96// "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
97//
a4e98857 98// Note, that one can set a title for new digits branch, and repeat digitization with
bca3b32a 99// another set of parameters.
100//
a4e98857 101// Use case:
990119d6 102// root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
103// root[1] d->ExecuteTask()
104// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
105// //Digitizes SDigitis in all events found in file galice.root
bca3b32a 106//
8cb3533f 107// root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;
108// // Will read sdigits from galice1.root
109// root[3] d1->MixWith("galice2.root")
990119d6 110// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
a4e98857 111// // Reads another set of sdigits from galice2.root
8cb3533f 112// root[3] d1->MixWith("galice3.root")
a4e98857 113// // Reads another set of sdigits from galice3.root
8cb3533f 114// root[4] d->ExecuteTask("deb timing")
115// // Reads SDigits from files galice1.root, galice2.root ....
116// // mixes them and stores produced Digits in file galice1.root
117// // deb - prints number of produced digits
118// // deb all - prints list of produced digits
119// // timing - prints time used for digitization
990119d6 120//
990119d6 121
122// --- ROOT system ---
990119d6 123#include "TTree.h"
124#include "TSystem.h"
8cb3533f 125#include "TBenchmark.h"
e957fea8 126#include "TRandom.h"
f898e0f3 127#include "TMath.h"
ba54256b 128
990119d6 129// --- Standard library ---
130
131// --- AliRoot header files ---
e250de8f 132#include <TGeoManager.h>
351dd634 133#include "AliLog.h"
3f81a70b 134#include "AliRunDigitizer.h"
990119d6 135#include "AliPHOSDigit.h"
990119d6 136#include "AliPHOSDigitizer.h"
8cb3533f 137#include "AliPHOSGeometry.h"
7437a0f7 138#include "AliPHOSTick.h"
f898e0f3 139#include "AliPHOSSimParam.h"
140#include "AliPHOSCalibData.h"
141#include "AliRunLoader.h"
142#include "AliPHOSLoader.h"
990119d6 143
144ClassImp(AliPHOSDigitizer)
145
146
147//____________________________________________________________________________
3663622c 148AliPHOSDigitizer::AliPHOSDigitizer() :
149 AliDigitizer("",""),
150 fDefaultInit(kTRUE),
151 fDigitsInRun(0),
152 fInit(kFALSE),
153 fInput(0),
154 fInputFileNames(0x0),
155 fEventNames(0x0),
156 fEmcCrystals(0),
3663622c 157 fEventFolderName(""),
158 fFirstEvent(0),
ddd1a39c 159 fLastEvent(0),
f898e0f3 160 fcdb(0x0),
ddd1a39c 161 fEventCounter(0)
990119d6 162{
163 // ctor
8d0f3f77 164 InitParameters() ;
fbf811ec 165 fManager = 0 ; // We work in the standalong mode
88cb7938 166}
3f81a70b 167
168//____________________________________________________________________________
e191bb57 169AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
170 TString eventFolderName):
3663622c 171 AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
172 fDefaultInit(kFALSE),
173 fDigitsInRun(0),
174 fInit(kFALSE),
175 fInput(0),
176 fInputFileNames(0x0),
177 fEventNames(0x0),
178 fEmcCrystals(0),
3663622c 179 fEventFolderName(eventFolderName),
180 fFirstEvent(0),
ddd1a39c 181 fLastEvent(0),
f898e0f3 182 fcdb(0x0),
ddd1a39c 183 fEventCounter(0)
3f81a70b 184{
185 // ctor
8d0f3f77 186 InitParameters() ;
3f81a70b 187 Init() ;
92f521a9 188 fDefaultInit = kFALSE ;
88cb7938 189 fManager = 0 ; // We work in the standalong mode
f898e0f3 190 fcdb = new AliPHOSCalibData(-1);
88cb7938 191}
192
193//____________________________________________________________________________
3663622c 194AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) :
195 AliDigitizer(d),
196 fDefaultInit(d.fDefaultInit),
197 fDigitsInRun(d.fDigitsInRun),
198 fInit(d.fInit),
199 fInput(d.fInput),
200 fInputFileNames(0x0),//?
201 fEventNames(0x0),//?
202 fEmcCrystals(d.fEmcCrystals),
3663622c 203 fEventFolderName(d.fEventFolderName),
204 fFirstEvent(d.fFirstEvent),
ddd1a39c 205 fLastEvent(d.fLastEvent),
f898e0f3 206 fcdb (0x0),
ddd1a39c 207 fEventCounter(0)
88cb7938 208{
209 // copyy ctor
88cb7938 210 SetName(d.GetName()) ;
211 SetTitle(d.GetTitle()) ;
f898e0f3 212 fcdb = new AliPHOSCalibData(-1);
990119d6 213}
214
990119d6 215//____________________________________________________________________________
3663622c 216AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
217 AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
218 fDefaultInit(kFALSE),
219 fDigitsInRun(0),
220 fInit(kFALSE),
221 fInput(0),
222 fInputFileNames(0x0),
223 fEventNames(0x0),
224 fEmcCrystals(0),
3663622c 225 fEventFolderName(fManager->GetInputFolderName(0)),
226 fFirstEvent(0),
ddd1a39c 227 fLastEvent(0),
f898e0f3 228 fcdb (0x0),
ddd1a39c 229 fEventCounter(0)
230
990119d6 231{
45fa49ca 232 // ctor Init() is called by RunDigitizer
88cb7938 233 fManager = rd ;
88cb7938 234 SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
b22e4735 235 InitParameters() ;
fbf811ec 236 fDefaultInit = kFALSE ;
f898e0f3 237 fcdb = new AliPHOSCalibData(-1);
990119d6 238}
239
240//____________________________________________________________________________
241 AliPHOSDigitizer::~AliPHOSDigitizer()
242{
243 // dtor
f898e0f3 244 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
245 if(rl){
246 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
247 phosLoader->CleanDigitizer() ;
248 }
249
88cb7938 250 delete [] fInputFileNames ;
251 delete [] fEventNames ;
ddd1a39c 252
f898e0f3 253 if(fcdb){ delete fcdb ; fcdb=0;}
ddd1a39c 254
990119d6 255}
256
257//____________________________________________________________________________
fc7e2f43 258void AliPHOSDigitizer::Digitize(Int_t event)
a4e98857 259{
260
261 // Makes the digitization of the collected summable digits.
262 // It first creates the array of all PHOS modules
dc986a1d 263 // filled with noise (different for EMC, and CPV) and
a4e98857 264 // then adds contributions from SDigits.
265 // This design avoids scanning over the list of digits to add
266 // contribution to new SDigits only.
990119d6 267
f898e0f3 268 Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
269
270 //First stream
271 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
272 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
273
274 Int_t readEvent = event ;
45fa49ca 275 if (fManager)
f898e0f3 276 readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
d6e8d7d3 277 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
f898e0f3 278 readEvent, GetTitle(), fEventFolderName.Data())) ;
279 rl->GetEvent(readEvent) ;
280 phosLoader->CleanSDigits() ;
281 phosLoader->LoadSDigits("READ") ;
282
283 //Prepare Output
284 TClonesArray * digits = phosLoader->Digits() ;
285 if( !digits ) {
286 phosLoader->MakeDigitsArray() ;
287 digits = phosLoader->Digits() ;
288 }
289
7b7c1533 290 digits->Clear() ;
990119d6 291
f898e0f3 292 //
293 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
990119d6 294 //Making digits with noise, first EMC
e250de8f 295 //Check which PHOS modules are present
296 Bool_t isPresent[5] ;
297 TString volpath ;
298 Int_t nmod=0 ;
299 for(Int_t i=0; i<5; i++){
300 volpath = "/ALIC_1/PHOS_";
301 volpath += i+1;
302 if (gGeoManager->CheckPath(volpath.Data())) {
303 isPresent[i]=1 ;
304 nmod++ ;
305 }
306 else{
307 isPresent[i]=0 ;
308 }
309 }
310
311 Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
990119d6 312
313 Int_t nCPV ;
990119d6 314 Int_t absID ;
990119d6 315
e250de8f 316 //check if CPV exists
317 Bool_t isCPVpresent=0 ;
318 for(Int_t i=1; i<=5 && !isCPVpresent; i++){
319 volpath = "/ALIC_1/PHOS_";
320 volpath += i;
321 volpath += "/PCPV_1";
322 if (gGeoManager->CheckPath(volpath.Data()))
323 isCPVpresent=1 ;
324 }
88cb7938 325
e250de8f 326 if(isCPVpresent){
327 nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
328 }
329 else{
330 nCPV = nEMC ;
331 }
332
9688c1dd 333 digits->Expand(nCPV) ;
8cb3533f 334
88cb7938 335 //take all the inputs to add together and load the SDigits
336 TObjArray * sdigArray = new TObjArray(fInput) ;
f898e0f3 337 sdigArray->AddAt(phosLoader->SDigits(), 0) ;
338
339 for(Int_t i = 1 ; i < fInput ; i++){
88cb7938 340 TString tempo(fEventNames[i]) ;
341 tempo += i ;
f898e0f3 342 AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
343 if(!rl2){
344 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
345 if (!rl2) {
346 Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
347 return ;
348 }
349 rl2->LoadHeader();
350 }
351 AliPHOSLoader * phosLoader2 = dynamic_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
352
353 if(fManager){
354 readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
355 }
356 TClonesArray * digs ;
357 if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
358 AliInfo(Form("Adding event %d from input stream %d %s %s",
359 readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
360 rl2->GetEvent(readEvent) ;
361 phosLoader2->CleanDigits() ;
362 phosLoader2->LoadDigits("READ") ;
363 digs = phosLoader2->Digits() ;
364 toMakeNoise=0 ; //Do not add noise, it is already in stream
365 }
366 else{
367 AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
368 readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
369 rl2->GetEvent(readEvent) ;
370 phosLoader2->CleanSDigits() ;
371 phosLoader2->LoadSDigits("READ") ;
372 digs = phosLoader2->SDigits() ;
373 }
374 sdigArray->AddAt(digs, i) ;
38bb0fd5 375 }
9688c1dd 376
27a73a5d 377 //Find the first crystal with signal
9688c1dd 378 Int_t nextSig = 200000 ;
88cb7938 379 TClonesArray * sdigits ;
f898e0f3 380 for(Int_t i = 0 ; i < fInput ; i++){
88cb7938 381 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
a6eedfad 382 if ( !sdigits->GetEntriesFast() )
7a9d98f9 383 continue ;
88cb7938 384 Int_t curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
385 if(curNext < nextSig)
386 nextSig = curNext ;
9688c1dd 387 }
88cb7938 388
389 TArrayI index(fInput) ;
9688c1dd 390 index.Reset() ; //Set all indexes to zero
88cb7938 391
9688c1dd 392 AliPHOSDigit * digit ;
393 AliPHOSDigit * curSDigit ;
88cb7938 394
f898e0f3 395// TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
88cb7938 396
9688c1dd 397 //Put Noise contribution
f898e0f3 398 Double_t apdNoise = 0. ;
399 if(toMakeNoise)
400 apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ;
401
e250de8f 402 Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
403 Int_t idigit= 0;
404 for(Int_t imod=0; imod<5; imod++){
405 if(!isPresent[imod])
406 continue ;
407 Int_t firstAbsId=imod*emcpermod+1 ;
408 Int_t lastAbsId =(imod+1)*emcpermod ;
409 for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
410 Float_t noise = gRandom->Gaus(0.,apdNoise) ;
411 new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
412 //look if we have to add signal?
413 digit = dynamic_cast<AliPHOSDigit *>(digits->At(idigit)) ;
414 idigit++ ;
88cb7938 415
e250de8f 416 if(absID==nextSig){
9688c1dd 417 //Add SDigits from all inputs
f898e0f3 418// ticks->Clear() ;
419// Int_t contrib = 0 ;
420
421//New Timing model is necessary
422// Float_t a = digit->GetEnergy() ;
423// Float_t b = TMath::Abs( a / fTimeSignalLength) ;
424// //Mark the beginning of the signal
425// new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
426// //Mark the end of the signal
427// new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
428
429// Calculate time as time of the largest digit
e250de8f 430 Float_t time = digit->GetTime() ;
431 Float_t eTime= digit->GetEnergy() ;
88cb7938 432
e250de8f 433 //loop over inputs
434 for(Int_t i = 0 ; i < fInput ; i++){
737c795f 435 if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
e250de8f 436 curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
737c795f 437 if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
438 curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
439 }
440 }
9688c1dd 441 else
e250de8f 442 curSDigit = 0 ;
443 //May be several digits will contribute from the same input
444 while(curSDigit && curSDigit->GetId() == absID){
445 //Shift primary to separate primaries belonging different inputs
446 Int_t primaryoffset ;
447 if(fManager)
448 primaryoffset = fManager->GetMask(i) ;
449 else
450 primaryoffset = 10000000*i ;
451 curSDigit->ShiftPrimary(primaryoffset) ;
7437a0f7 452
f898e0f3 453//New Timing model is necessary
454// a = curSDigit->GetEnergy() ;
455// b = a /fTimeSignalLength ;
456// new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
457// new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
458 if(curSDigit->GetEnergy()>eTime){
459 eTime=curSDigit->GetEnergy() ;
460 time=curSDigit->GetTime() ;
461 }
e250de8f 462 *digit += *curSDigit ; //add energies
463
464 index[i]++ ;
465 if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
466 curSDigit = dynamic_cast<AliPHOSDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
467 else
468 curSDigit = 0 ;
469 }
470 }
88cb7938 471
e250de8f 472 //calculate and set time
f898e0f3 473//New Timing model is necessary
474// Float_t time = FrontEdgeTime(ticks) ;
e250de8f 475 digit->SetTime(time) ;
88cb7938 476
e250de8f 477 //Find next signal module
478 nextSig = 200000 ;
479 for(Int_t i = 0 ; i < fInput ; i++){
480 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
481 Int_t curNext = nextSig ;
482 if(sdigits->GetEntriesFast() > index[i] ){
483 curNext = dynamic_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
484 }
485 if(curNext < nextSig) nextSig = curNext ;
486 }
7b7c1533 487 }
488 }
990119d6 489 }
f898e0f3 490
491 //distretize energy if necessary
492 if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
493 Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
494 for(Int_t i = 0 ; i < nEMC ; i++){
495 digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
496 digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
497 }
498 }
737c795f 499 //Apply decalibration if necessary
500 for(Int_t i = 0 ; i < nEMC ; i++){
501 digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
502 Decalibrate(digit) ;
167d2780 503 }
504
3f81a70b 505
f898e0f3 506// ticks->Delete() ;
507// delete ticks ;
88cb7938 508
9688c1dd 509 //Now CPV digits (different noise and no timing)
e250de8f 510 Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
511 Int_t nEMCtotal=emcpermod*5 ;
f898e0f3 512 Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
e250de8f 513 if(isCPVpresent){ //CPV is present in current geometry
514 for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
515 if(!isPresent[imod])
516 continue ;
517 Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
518 Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
519 for(absID = firstAbsId; absID <= lastAbsId; absID++){
520 Float_t noise = gRandom->Gaus(0., cpvNoise) ;
521 new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
522 idigit++ ;
523 //look if we have to add signal?
524 if(absID==nextSig){
525 digit = dynamic_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
526 //Add SDigits from all inputs
527 for(Int_t i = 0 ; i < fInput ; i++){
528 if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
529 curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
530 else
531 curSDigit = 0 ;
532
533 //May be several digits will contribute from the same input
534 while(curSDigit && curSDigit->GetId() == absID){
535 //Shift primary to separate primaries belonging different inputs
536 Int_t primaryoffset ;
537 if(fManager)
538 primaryoffset = fManager->GetMask(i) ;
539 else
540 primaryoffset = 10000000*i ;
541 curSDigit->ShiftPrimary(primaryoffset) ;
542
543 //add energies
544 *digit += *curSDigit ;
545 index[i]++ ;
546 if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
547 curSDigit = dynamic_cast<AliPHOSDigit*>( dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
548 else
549 curSDigit = 0 ;
550 }
551 }
a6eedfad 552
e250de8f 553 //Find next signal module
554 nextSig = 200000 ;
555 for(Int_t i = 0 ; i < fInput ; i++){
556 sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
557 Int_t curNext = nextSig ;
558 if(sdigits->GetEntriesFast() > index[i] )
559 curNext = dynamic_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
560 if(curNext < nextSig) nextSig = curNext ;
561 }
9688c1dd 562
e250de8f 563 }
564 }
9688c1dd 565 }
566 }
88cb7938 567
568 delete sdigArray ; //We should not delete its contents
9688c1dd 569
1dfe0f3c 570 Int_t relId[4];
571
572 //set amplitudes in bad channels to zero
f898e0f3 573
65fd5dc1 574 for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
1dfe0f3c 575 digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
f898e0f3 576 geom->AbsToRelNumbering(digit->GetId(),relId);
1dfe0f3c 577 if(relId[1] == 0) // Emc
f898e0f3 578 if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
1dfe0f3c 579 }
580
990119d6 581 //remove digits below thresholds
f898e0f3 582 Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
583 for(Int_t i = 0 ; i < nEMC ; i++){
548f0134 584 digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
a28333c5 585
f898e0f3 586 if(digit->GetEnergy() < emcThreshold){
a6eedfad 587 digits->RemoveAt(i) ;
f898e0f3 588 continue ;
589 }
590 Float_t tres = TimeResolution(digit->GetEnergy()) ;
591 digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
aaf8a71c 592 }
593
f898e0f3 594 Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
595 for(Int_t i = nEMC; i < nCPV ; i++){
596 if( dynamic_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
a6eedfad 597 digits->RemoveAt(i) ;
f898e0f3 598 }
9688c1dd 599
7b7c1533 600 digits->Compress() ;
7b7c1533 601 Int_t ndigits = digits->GetEntriesFast() ;
a91822a8 602
3758d9fc 603 //Set indexes in list of digits and make true digitization of the energy
f898e0f3 604 for (Int_t i = 0 ; i < ndigits ; i++) {
548f0134 605 digit = dynamic_cast<AliPHOSDigit*>( digits->At(i) ) ;
990119d6 606 digit->SetIndexInList(i) ;
27a73a5d 607 if(digit->GetId() > fEmcCrystals){ //digitize CPV only
877695e7 608 digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
27a73a5d 609 }
990119d6 610 }
f898e0f3 611
7b7c1533 612}
737c795f 613//____________________________________________________________________________
614Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
615 //Apply calibration
616 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
548f0134 617
737c795f 618 //Determine rel.position of the cell absolute ID
619 Int_t relId[4];
620 geom->AbsToRelNumbering(absId,relId);
621 Int_t module=relId[0];
622 Int_t row =relId[2];
623 Int_t column=relId[3];
624 if(relId[1]==0){ //This Is EMC
625 Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
626 return amp*calibration ;
627 }
628 return 0 ;
629}
3758d9fc 630//____________________________________________________________________________
737c795f 631void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
3758d9fc 632{
877695e7 633 // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
634
f898e0f3 635 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
877695e7 636
637 //Determine rel.position of the cell absolute ID
638 Int_t relId[4];
f898e0f3 639 geom->AbsToRelNumbering(digit->GetId(),relId);
877695e7 640 Int_t module=relId[0];
641 Int_t row =relId[2];
642 Int_t column=relId[3];
737c795f 643 if(relId[1]==0){ //This Is EMC
644 Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
645 Float_t energy = digit->GetEnergy()/calibration;
646 digit->SetEnergy(energy); //Now digit measures E in ADC counts
647 }
877695e7 648}
649//____________________________________________________________________________
650Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
651{
652 // Returns digitized value of the CPV charge in a pad absId
0bc3b8ed 653
f898e0f3 654 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
fc6706cb 655
a8ec0771 656 //Determine rel.position of the cell absId
657 Int_t relId[4];
f898e0f3 658 geom->AbsToRelNumbering(absId,relId);
a8ec0771 659 Int_t module=relId[0];
27a73a5d 660 Int_t row =relId[2];
a8ec0771 661 Int_t column=relId[3];
662
877695e7 663 Int_t channel = 0;
a8ec0771 664
877695e7 665 if(absId > fEmcCrystals){ //digitize CPV only
a8ec0771 666
667 //reading calibration data for cell absId.
f898e0f3 668 Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
669 Float_t adcChanelCpv = fcdb->GetADCchannelCpv( module,column,row);
a8ec0771 670
f898e0f3 671 channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
672 Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
673 if(channel > nMax ) channel = nMax ;
3758d9fc 674 }
877695e7 675 return channel ;
3758d9fc 676}
548f0134 677
7b7c1533 678//____________________________________________________________________________
679void AliPHOSDigitizer::Exec(Option_t *option)
680{
212d1c0f 681 // Steering method to process digitization for events
682 // in the range from fFirstEvent to fLastEvent.
683 // This range is optionally set by SetEventRange().
45fa49ca 684 // if fLastEvent=-1, then process events until the end.
685 // by default fLastEvent = fFirstEvent (process only one event)
88cb7938 686
687 if (!fInit) { // to prevent overwrite existing file
351dd634 688 AliError(Form("Give a version name different from %s",
689 fEventFolderName.Data() )) ;
88cb7938 690 return ;
691 }
990119d6 692
7b7c1533 693 if (strstr(option,"print")) {
88cb7938 694 Print();
7b7c1533 695 return ;
8cb3533f 696 }
990119d6 697
8661738e 698 if(strstr(option,"tim"))
699 gBenchmark->Start("PHOSDigitizer");
3f81a70b 700
f898e0f3 701 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
702 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
88cb7938 703
5b4d2c50 704 // Post Digitizer to the white board
f898e0f3 705 phosLoader->PostDigitizer(this) ;
5b4d2c50 706
212d1c0f 707 if (fLastEvent == -1)
f898e0f3 708 fLastEvent = rl->GetNumberOfEvents() - 1 ;
396a348e 709 else if (fManager)
710 fLastEvent = fFirstEvent ;
45fa49ca 711
212d1c0f 712 Int_t nEvents = fLastEvent - fFirstEvent + 1;
713
7b7c1533 714 Int_t ievent ;
88cb7938 715
212d1c0f 716 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
ddd1a39c 717 fEventCounter++ ;
b22e4735 718
7b7c1533 719 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
88cb7938 720
90cceaf6 721 WriteDigits() ;
88cb7938 722
01a599c9 723 if(strstr(option,"deb"))
724 PrintDigits(option);
94de8339 725
726 //increment the total number of Digits per run
f898e0f3 727 fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
88cb7938 728 }
ddd1a39c 729
f898e0f3 730 phosLoader->CleanDigitizer();
5b4d2c50 731
8cb3533f 732 if(strstr(option,"tim")){
733 gBenchmark->Stop("PHOSDigitizer");
21cd0c07 734 TString message ;
735 message = " took %f seconds for Digitizing %f seconds per event\n" ;
351dd634 736 AliInfo(Form( message.Data(),
21cd0c07 737 gBenchmark->GetCpuTime("PHOSDigitizer"),
351dd634 738 gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
21cd0c07 739 }
990119d6 740}
7b7c1533 741//____________________________________________________________________________
f898e0f3 742Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
743 //calculate TOF resolution using beam-test resutls
744 Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
745 Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
746 return TMath::Sqrt(a*a+b*b/e) ;
9688c1dd 747}
8d0f3f77 748
f898e0f3 749////____________________________________________________________________________
750//Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
751//{
752// // Returns the shortest time among all time ticks
753//
754// ticks->Sort() ; //Sort in accordance with times of ticks
755// TIter it(ticks) ;
756// AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
757// Float_t time = ctick->CrossingTime(fTimeThreshold) ;
758//
759// AliPHOSTick * t ;
760// while((t=(AliPHOSTick*) it.Next())){
761// if(t->GetTime() < time) //This tick starts before crossing
762// *ctick+=*t ;
763// else
764// return time ;
765//
766// time = ctick->CrossingTime(fTimeThreshold) ;
767// }
768// return time ;
769//}
770
9688c1dd 771//____________________________________________________________________________
3f81a70b 772Bool_t AliPHOSDigitizer::Init()
a4e98857 773{
fbf811ec 774 // Makes all memory allocations
88cb7938 775 fInit = kTRUE ;
f898e0f3 776
bfae5a5d 777 AliPHOSGeometry *geom;
778 if (!(geom = AliPHOSGeometry::GetInstance()))
779 geom = AliPHOSGeometry::GetInstance("IHEP","");
780// const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
b22e4735 781
88cb7938 782 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
8d0f3f77 783
45fa49ca 784 fFirstEvent = 0 ;
785 fLastEvent = fFirstEvent ;
88cb7938 786 if (fManager)
787 fInput = fManager->GetNinputs() ;
788 else
789 fInput = 1 ;
790
791 fInputFileNames = new TString[fInput] ;
792 fEventNames = new TString[fInput] ;
793 fInputFileNames[0] = GetTitle() ;
794 fEventNames[0] = fEventFolderName.Data() ;
795 Int_t index ;
796 for (index = 1 ; index < fInput ; index++) {
797 fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
798 TString tempo = fManager->GetInputFolderName(index) ;
45fa49ca 799 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
8d0f3f77 800 }
88cb7938 801
802 //to prevent cleaning of this object while GetEvent is called
f898e0f3 803 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
804 if(!rl){
805 rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ;
806 }
807 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
808 phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
88cb7938 809
810 return fInit ;
8d0f3f77 811}
812
813//____________________________________________________________________________
814void AliPHOSDigitizer::InitParameters()
815{
45fa49ca 816 // Set initial parameters Digitizer
0bc3b8ed 817
3758d9fc 818 fDigitsInRun = 0 ;
212d1c0f 819 SetEventRange(0,-1) ;
8cb3533f 820
990119d6 821}
7b7c1533 822
990119d6 823//__________________________________________________________________
702ab87e 824void AliPHOSDigitizer::Print(const Option_t *)const
21cd0c07 825{
dd5c4038 826 // Print Digitizer's parameters
351dd634 827 AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
88cb7938 828 if( strcmp(fEventFolderName.Data(), "") != 0 ){
829 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
830
831 Int_t nStreams ;
832 if (fManager)
833 nStreams = GetNInputStreams() ;
834 else
835 nStreams = fInput ;
836
837 Int_t index = 0 ;
838 for (index = 0 ; index < nStreams ; index++) {
839 TString tempo(fEventNames[index]) ;
840 tempo += index ;
f898e0f3 841 printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
88cb7938 842 }
88cb7938 843
f898e0f3 844 // printf("\nWith following parameters:\n") ;
845 // printf(" Electronics noise in EMC (fPinNoise) = %f GeV\n", fPinNoise ) ;
846 // printf(" Threshold in EMC (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
847 // printf(" Noise in CPV (fCPVNoise) = %f aux units\n", fCPVNoise ) ;
848 // printf(" Threshold in CPV (fCPVDigitThreshold) = %f aux units\n",fCPVDigitThreshold ) ;
88cb7938 849 printf(" ---------------------------------------------------\n") ;
990119d6 850 }
8cb3533f 851 else
351dd634 852 AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
8cb3533f 853
854}
88cb7938 855
8cb3533f 856//__________________________________________________________________
21cd0c07 857 void AliPHOSDigitizer::PrintDigits(Option_t * option)
858{
3bf72d32 859 // Print a table of digits
21cd0c07 860
f898e0f3 861
862 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
863 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
864 TClonesArray * digits = phosLoader->Digits() ;
865 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
3bf72d32 866
351dd634 867 AliInfo(Form("%d", digits->GetEntriesFast())) ;
88cb7938 868 printf("\nevent %d", gAlice->GetEvNumber()) ;
869 printf("\n Number of entries in Digits list %d", digits->GetEntriesFast() ) ;
870
11f9c5ff 871
3bf72d32 872 if(strstr(option,"all")||strstr(option,"EMC")){
8cb3533f 873 //loop over digits
874 AliPHOSDigit * digit;
88cb7938 875 printf("\nEMC digits (with primaries):\n") ;
876 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
f898e0f3 877 Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
8cb3533f 878 Int_t index ;
a6eedfad 879 for (index = 0 ; (index < digits->GetEntriesFast()) &&
88cb7938 880 (dynamic_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
7b7c1533 881 digit = (AliPHOSDigit * ) digits->At(index) ;
21cd0c07 882 if(digit->GetNprimary() == 0)
883 continue;
27a73a5d 884// printf("%6d %8d %6.5e %4d %2d :",
885// digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ; // YVK
886 printf("%6d %.4f %6.5e %4d %2d :",
887 digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
8cb3533f 888 Int_t iprimary;
21cd0c07 889 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
88cb7938 890 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
21cd0c07 891 }
81d5f731 892 printf("\n") ;
21cd0c07 893 }
9688c1dd 894 }
3bf72d32 895
9688c1dd 896 if(strstr(option,"all")||strstr(option,"CPV")){
8cb3533f 897
9688c1dd 898 //loop over CPV digits
899 AliPHOSDigit * digit;
88cb7938 900 printf("\nCPV digits (with primaries):\n") ;
901 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
f898e0f3 902 Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
9688c1dd 903 Int_t index ;
a6eedfad 904 for (index = 0 ; index < digits->GetEntriesFast(); index++) {
9688c1dd 905 digit = (AliPHOSDigit * ) digits->At(index) ;
906 if(digit->GetId() > maxEmc){
81d5f731 907 printf("%6d %8d %4d %2d :",
11f9c5ff 908 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
9688c1dd 909 Int_t iprimary;
21cd0c07 910 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
88cb7938 911 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
21cd0c07 912 }
81d5f731 913 printf("\n") ;
21cd0c07 914 }
9688c1dd 915 }
8cb3533f 916 }
88cb7938 917
8cb3533f 918}
7b7c1533 919
9688c1dd 920//__________________________________________________________________
0bc3b8ed 921Float_t AliPHOSDigitizer::TimeOfNoise(void) const
9688c1dd 922{ // Calculates the time signal generated by noise
26a2ef9d 923 //PH Info("TimeOfNoise", "Change me") ;
04f0bda3 924 return gRandom->Rndm() * 1.28E-5;
9688c1dd 925}
7b7c1533 926
88cb7938 927//__________________________________________________________________
928void AliPHOSDigitizer::Unload()
929{
930
931 Int_t i ;
932 for(i = 1 ; i < fInput ; i++){
933 TString tempo(fEventNames[i]) ;
934 tempo += i ;
f898e0f3 935 AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
936 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
937 phosLoader->UnloadSDigits() ;
88cb7938 938 }
939
f898e0f3 940 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
941 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
942 phosLoader->UnloadDigits() ;
88cb7938 943}
944
7b7c1533 945//____________________________________________________________________________
90cceaf6 946void AliPHOSDigitizer::WriteDigits()
7b7c1533 947{
948
949 // Makes TreeD in the output file.
950 // Check if branch already exists:
951 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
952 // already existing branches.
953 // else creates branch with Digits, named "PHOS", title "...",
954 // and branch "AliPHOSDigitizer", with the same title to keep all the parameters
955 // and names of files, from which digits are made.
956
f898e0f3 957 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
958 AliPHOSLoader * phosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
959
960 const TClonesArray * digits = phosLoader->Digits() ;
961 TTree * treeD = phosLoader->TreeD();
962 if(!treeD){
963 phosLoader->MakeTree("D");
964 treeD = phosLoader->TreeD();
965 }
966
7b7c1533 967 // -- create Digits branch
968 Int_t bufferSize = 32000 ;
2524c56f 969 TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
88cb7938 970 digitsBranch->SetTitle(fEventFolderName);
971 digitsBranch->Fill() ;
fbf811ec 972
f898e0f3 973 phosLoader->WriteDigits("OVERWRITE");
974 phosLoader->WriteDigitizer("OVERWRITE");
88cb7938 975
976 Unload() ;
b3690abb 977
7b7c1533 978}