EMCAL/AliEMCALSDigitizer.cxx

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