1 /*************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 in one tower of the EMCAL
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.
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
44 //*-- Author : Sahal Yacoob (LBL)
45 // based on : AliPHOSSDigitzer
47 // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
48 // of new IO (à la PHOS)
49 //////////////////////////////////////////////////////////////////////////////
51 // --- ROOT system ---
58 #include "TBenchmark.h"
59 #include "TGeometry.h"
61 // --- Standard library ---
63 // --- AliRoot header files ---
65 #include "AliHeader.h"
66 #include "AliEMCALDigit.h"
67 #include "AliEMCALGeometry.h"
68 #include "AliEMCALGetter.h"
69 #include "AliEMCALHit.h"
70 #include "AliEMCALSDigitizer.h"
72 ClassImp(AliEMCALSDigitizer)
74 //____________________________________________________________________________
75 AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
81 fDefaultInit = kTRUE ;
84 //____________________________________________________________________________
85 AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle, const Bool_t toSplit):
86 TTask(sDigitsTitle, headerFile)
92 fDefaultInit = kFALSE ;
95 //____________________________________________________________________________
96 AliEMCALSDigitizer::~AliEMCALSDigitizer()
102 //____________________________________________________________________________
103 void AliEMCALSDigitizer::Init(){
104 // Initialization: open root-file, allocate arrays for hits and sdigits,
105 // attach task SDigitizer to the list of EMCAL tasks
107 // Initialization can not be done in the default constructor
108 //============================================================= YS
109 // The initialisation is now done by the getter
111 if( strcmp(GetTitle(), "") == 0 )
112 SetTitle("galice.root") ;
114 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), fToSplit) ;
116 Error("Init", "Could not obtain the Getter object !" ) ;
120 gime->PostSDigits( GetName(), GetTitle() ) ;
125 // construct the name of the file as /path/EMCAL.SDigits.root
126 // First - extract full path if necessary
127 TString sDigitsFileName(GetTitle()) ;
128 Ssiz_t islash = sDigitsFileName.Last('/') ;
129 if(islash<sDigitsFileName.Length())
130 sDigitsFileName.Remove(islash+1,sDigitsFileName.Length()) ;
134 // Next - append the file name
135 sDigitsFileName+="EMCAL.SDigits." ;
136 if((strcmp(GetName(),"Default")!=0)&&(strcmp(GetName(),"")!=0)){
137 sDigitsFileName+=GetName() ;
138 sDigitsFileName+="." ;
140 sDigitsFileName+="root" ;
142 // Finally - check if the file already opened or open the file
143 fSplitFile = static_cast<TFile*>(gROOT->GetFile(sDigitsFileName.Data()));
145 fSplitFile = TFile::Open(sDigitsFileName.Data(),"update") ;
148 TString sdname(GetName() );
150 sdname.Append(GetTitle() ) ;
152 gime->PostSDigitizer(this) ;
156 //____________________________________________________________________________
157 void AliEMCALSDigitizer::InitParameters()
162 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
163 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
164 if (geom->GetSampling() == 0.) {
165 Error("InitParameters", "Sampling factor not set !") ;
169 Info("InitParameters", "Sampling factor set to %f\n", geom->GetSampling()) ;
171 // this threshold corresponds approximately to 100 MeV
172 fECPrimThreshold = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNECLayers() ;
173 fPREPrimThreshold = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNPRLayers() ;
174 fHCPrimThreshold = fECPrimThreshold/5. ; // 5 is totally arbitrary
178 //____________________________________________________________________________
179 void AliEMCALSDigitizer::Exec(Option_t *option)
181 // Collects all hits in the section (PRE/ECAL/HCAL) of the same tower into digit
183 if( strcmp(GetName(), "") == 0 )
186 if (strstr(option, "print") ) {
191 if(strstr(option,"tim"))
192 gBenchmark->Start("EMCALSDigitizer");
194 //Check, if this branch already exits
195 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
196 if(gime->BranchExists("SDigits") )
199 TString sdname(GetName()) ;
200 sdname.Remove(sdname.Index(GetTitle())-1) ;
202 Int_t nevents = gime->MaxEvent() ;
204 for(ievent = 0; ievent < nevents; ievent++){
205 gime->Event(ievent,"H") ;
206 const TClonesArray * hits = gime->Hits() ;
207 TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
211 //Now make SDigits from hits, for EMCAL it is the same, so just copy
212 Int_t nPrim = static_cast<Int_t>((gAlice->TreeH())->GetEntries()) ;
213 // Attention nPrim is the number of primaries tracked by Geant
214 // and this number could be different to the number of Primaries in TreeK;
216 for ( iprim = 0 ; iprim < nPrim ; iprim++ ) {
217 //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
220 for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
221 AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
222 AliEMCALDigit * curSDigit = 0 ;
223 AliEMCALDigit * sdigit = 0 ;
224 Bool_t newsdigit = kTRUE;
226 // Assign primary number only if deposited energy is significant
228 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
231 Info("Exec", "id = %d energy = %f thresholdPRE = %f thresholdEC = %f thresholdHC = %f \n",
232 hit->GetId(), hit->GetEnergy(), fPREPrimThreshold, fECPrimThreshold, fHCPrimThreshold) ;
234 if( geom->IsInPRE(hit->GetId()) )
235 if( hit->GetEnergy() > fPREPrimThreshold )
236 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
237 hit->GetIparent(), hit->GetId(),
238 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
240 curSDigit = new AliEMCALDigit( -1 ,
243 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
244 else if( geom->IsInECAL(hit->GetId()) )
245 if( hit->GetEnergy() > fECPrimThreshold )
246 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
247 hit->GetIparent(), hit->GetId(),
248 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
250 curSDigit = new AliEMCALDigit( -1 ,
253 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
254 else if( geom->IsInHCAL(hit->GetId()) )
255 if( hit->GetEnergy() > fHCPrimThreshold )
257 curSDigit = new AliEMCALDigit( hit->GetPrimary(),
258 hit->GetIparent(), hit->GetId(),
259 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
261 curSDigit = new AliEMCALDigit( -1 ,
264 Digitize(hit->GetEnergy()), hit->GetTime() ) ;
267 for(check= 0; check < nSdigits ; check++) {
268 sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
269 if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL/HCAL/preshower tower ?
270 *sdigit = *sdigit + *curSDigit;
275 new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit);
279 } // loop over all hits (hit = deposited energy/layer/entering particle)
284 nSdigits = sdigits->GetEntriesFast() ;
285 fSDigitsInRun += nSdigits ;
286 sdigits->Expand(nSdigits) ;
288 Int_t NPrimarymax = -1 ;
290 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
291 AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
292 sdigit->SetIndexInList(i) ;
295 for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
296 if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
297 NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
302 if(gAlice->TreeS() == 0 || (fSplitFile)) //<--- To be checked: we should not create TreeS if it is already here
303 gAlice->MakeTree("S",fSplitFile);
308 //First list of sdigits
309 Int_t bufferSize = 32000 ;
310 TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
311 sdigitsBranch->SetTitle(sdname);
314 Int_t splitlevel = 0 ;
315 AliEMCALSDigitizer * sd = this ;
316 TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
317 &sd,bufferSize,splitlevel);
318 sdigitizerBranch->SetTitle(sdname);
320 sdigitsBranch->Fill() ;
321 sdigitizerBranch->Fill() ;
323 gAlice->TreeS()->AutoSave() ;
325 if(strstr(option,"deb"))
326 PrintSDigits(option) ;
329 if(strstr(option,"tim")){
330 gBenchmark->Stop("EMCALSDigitizer");
331 Info("Exec", "took %f seconds for SDigitizing %f seconds per event",
332 gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ;
336 //__________________________________________________________________
337 void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
339 // Setting title to branch SDigits
341 TString stitle(title) ;
343 // check if branch with title already exists
344 TBranch * sdigitsBranch =
345 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
346 TBranch * sdigitizerBranch =
347 static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
348 const char * sdigitsTitle = sdigitsBranch ->GetTitle() ;
349 const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
350 if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
351 Error("SetSDigitsBranch", "Cannot overwrite existing branch with title %s", title) ;
355 Info("SetSDigitsBranch", "Changing SDigits file from %s to %s", GetName(), title) ;
359 // Post to the WhiteBoard
360 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
361 gime->PostSDigits( title, GetTitle()) ;
365 //__________________________________________________________________
366 void AliEMCALSDigitizer::Print(Option_t* option)const
368 // Prints parameters of SDigitizer
370 TString message("\n") ;
371 message += "------------------- ";
372 message += GetName() ;
373 message += " -------------\n" ;
374 message += " Writing SDigitis to branch with title " ;
375 message += GetName() ;
376 message += "\n with digitization parameters A = " ;
378 message += "\n B = " ;
380 message += "\n Threshold for Primary assignment in PreShower = " ;
381 message += fPREPrimThreshold ;
382 message += "\n Threshold for Primary assignment in EC section= " ;
383 message += fECPrimThreshold ;
384 message += "\n Threshold for Primary assignment in HC section= " ;
385 message += fHCPrimThreshold ;
386 message += "\n---------------------------------------------------" ;
388 Info("Print", message.Data() ) ;
391 //__________________________________________________________________
392 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
395 // SDititizers are equal if their pedestal, slope and threshold are equal
397 if( (fA==sd.fA)&&(fB==sd.fB)&&
398 (fECPrimThreshold==sd.fECPrimThreshold) &&
399 (fHCPrimThreshold==sd.fHCPrimThreshold) &&
400 (fPREPrimThreshold==sd.fPREPrimThreshold))
406 //__________________________________________________________________
407 void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
408 //Prints list of digits produced at the current pass of AliEMCALDigitizer
410 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
411 TString sdname(GetName()) ;
412 sdname.Remove(sdname.Index(GetTitle())-1) ;
413 const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
415 TString message("\n") ;
416 message += "event " ;
417 message += gAlice->GetEvNumber() ;
418 message += "\n Number of entries in SDigits list " ;
419 message += sdigits->GetEntriesFast() ;
420 if(strstr(option,"all")||strstr(option,"EMC")){
423 AliEMCALDigit * digit;
424 message += "\n Id Amplitude Time Index Nprim: Primaries list \n" ;
426 char * tempo = new char[8192];
427 for (index = 0 ; index < sdigits->GetEntries() ; index++) {
428 digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
429 sprintf(tempo, "\n%6d %8d %6.5e %4d %2d :",
430 digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
434 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
435 sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
441 Info("PrintSDigits", message.Data() ) ;
444 //_______________________________________________________________________________________
445 // void AliEMCALSDigitizer::TestTowerID(void)
449 // Bool_t preshower = kFALSE;
450 // for (j = 0 ; j < 10 ; j++){ // loop over hit id
452 // for (i = 0 ; i <= 2 ; i++){ // loop over
453 // Int_t k = i*96*144+j*144+1;
454 // Info("TestTowerID", " Hit Index = %d %d TOWERID = %d", k, j*10, Layer2TowerID(k, preshower) ) ;
459 //____________________________________________________________________________
460 void AliEMCALSDigitizer::UseHitsFrom(const char * filename)