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 //////////////////////////////////////////////////////////////////////////////
  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"
  67 #include "AliEMCALGetter.h"
  68
  69 ClassImp(AliEMCALSDigitizer)
  70
  71
  72 //____________________________________________________________________________
  73   AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","")
  74 {
  75   // ctor
  76   fA = 0;
  77   fB = 10000000. ;
  78   fPrimThreshold = 0.01 ;
  79   fNevents = 0 ;
  80   fSDigits = 0 ;
  81   fHits = 0 ;
  82   fIsInitialized = kFALSE ;
  83
  84 }
  85
  86 //____________________________________________________________________________
  87 AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask(sDigitsTitle, headerFile)
  88 {
  89   // ctor
  90   fA = 0;
  91   fB = 10000000.;
  92   fPrimThreshold = 0.01 ;
  93   fNevents = 0 ;
  94   Init();
  95 }
  96
  97 //____________________________________________________________________________
  98 AliEMCALSDigitizer::~AliEMCALSDigitizer()
  99 {
 100   // dtor
 101   if(fSDigits)
 102     delete fSDigits ;
 103   if(fHits)
 104     delete fHits ;
 105 }
 106 //____________________________________________________________________________
 107 void AliEMCALSDigitizer::Init(){
 108
 109   // Initialization: open root-file, allocate arrays for hits and sdigits,
 110   // attach task SDigitizer to the list of PHOS tasks
 111   //
 112   // Initialization can not be done in the default constructor
 113   //============================================================= YS
 114   //  The initialisation is now done by the getter
 115
 116 if( strcmp(GetTitle(), "") == 0 )
 117     SetTitle("galice.root") ;
 118
 119   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), "update") ;
 120   if ( gime == 0 ) {
 121     cerr << "ERROR: AliEMCALSDigitizer::Init -> Could not obtain the Getter object !" << endl ;
 122     return ;
 123   }
 124
 125   gime->PostSDigits( GetName(), GetTitle() ) ;
 126
 127   TString sdname(GetName() );
 128   sdname.Append(":") ;
 129   sdname.Append(GetTitle() ) ;
 130   SetName(sdname) ;
 131   gime->PostSDigitizer(this) ;
 132
 133   // create a folder on the white board //YSAlice/WhiteBoard/SDigits/EMCAL/headerFile/sdigitsTitle
 134
 135  }
 136 //____________________________________________________________________________
 137 void AliEMCALSDigitizer::Exec(Option_t *option) {
 138
 139 // Collects all hits in the same active volume into digit
 140
 141   if( strcmp(GetName(), "") == 0 )
 142     Init() ;
 143
 144   if (strstr(option, "print") ) {
 145     Print("") ;
 146     return ;
 147   }
 148
 149   if(strstr(option,"tim"))
 150     gBenchmark->Start("EMCALSDigitizer");
 151   //Check, if this branch already exits
 152   gAlice->GetEvent(0) ;
 153   if(gAlice->TreeS() ) {
 154     TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
 155     TIter next(lob) ;
 156     TBranch * branch = 0 ;
 157     Bool_t emcalfound = kFALSE, sdigitizerfound = kFALSE ;
 158
 159     while ( (branch = (static_cast<TBranch*>(next()))) && (!emcalfound || !sdigitizerfound) ) {
 160       if ( (strcmp(branch->GetName(), "EMCAL")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
 161         emcalfound = kTRUE ;
 162
 163       else if ( (strcmp(branch->GetName(), "AliEMCALSDigitizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) )
 164         sdigitizerfound = kTRUE ;
 165     }
 166
 167     if ( emcalfound || sdigitizerfound ) {
 168       cerr << "WARNING: AliEMCALSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName()
 169            << " already exits" << endl ;
 170       return ;
 171     }
 172   }
 173   TString sdname(GetName()) ;
 174   sdname.Remove(sdname.Index(GetTitle())-1) ;
 175   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 176   const AliEMCALGeometry *geom = gime->EMCALGeometry();
 177   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
 178
 179   Int_t ievent ;
 180   for(ievent = 0; ievent < nevents; ievent++){
 181     gime->Event(ievent,"H") ;
 182     const TClonesArray * fHits = gime->Hits() ;
 183     TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
 184     sdigits->Clear();
 185     Int_t nSdigits = 0 ;
 186
 187   //Collects all hits in the same active volume into digit
 188
 189
 190
 191
 192     //Now made SDigits from hits, for EMCAL it is the same, so just copy
 193     //  Int_t itrack ;
 194     // for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
 195     //gime->Track(itrack);
 196       //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack
 197
 198
 199
 200
 201       Int_t i;
 202       for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
 203         AliEMCALHit * hit = (AliEMCALHit*)fHits->At(i) ;
 204         AliEMCALDigit  * curSDigit = 0  ;
 205         AliEMCALDigit * sdigit ;
 206         Bool_t newsdigit = kTRUE;
 207         // Assign primary number only if contribution is significant
 208         if( hit->GetEnergy() > fPrimThreshold)
 209           curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
 210                                           hit->GetIparent(),Layer2TowerID(hit->GetId(),kFALSE),
 211                                           Digitize( hit->GetEnergy() ) ,
 212                                           hit->GetTime()) ;
 213         else
 214           curSDigit =  new AliEMCALDigit( -1               ,
 215                                           -1               ,
 216                                           Layer2TowerID(hit->GetId(),kFALSE),
 217                                           Digitize( hit->GetEnergy() ) ,
 218                                           hit->GetTime() ) ;
 219         //cout << "Hit ID = " <<hit->GetId() << endl ;
 220         //cout << "ID for detector = " << curSDigit->GetId() << endl ;
 221         //  cout << hit->GetEnergy() << " - hit energy   -   Digit Energy - " << curSDigit->GetAmp() << endl;
 222         for(Int_t check= 0; check < nSdigits ; check++) {
 223           sdigit = (AliEMCALDigit *)sdigits->At(check);
 224           if( sdigit->GetId() == curSDigit->GetId())
 225             {// cout << "SDigit - Get Amp  " << sdigit->GetAmp() << endl ;
 226               *sdigit = *sdigit + *curSDigit ;
 227               newsdigit = kFALSE;
 228               //  cout << " and after addition " << sdigit->GetAmp() << endl ;
 229             }
 230         }
 231         if (newsdigit)
 232           { new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
 233           nSdigits++ ;
 234           //cout << "Detector nsdigits = " << nSdigits << endl ;
 235           }
 236         newsdigit = kTRUE;
 237
 238
 239         if( hit->GetEnergy() > fPrimThreshold)
 240           curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
 241                                           hit->GetIparent(),
 242                                           Layer2TowerID(hit->GetId(),kTRUE) ,
 243                                           Digitize( hit->GetEnergy() ),
 244                                           hit->GetTime( ) ) ;
 245         else
 246           curSDigit =  new AliEMCALDigit( -1               ,
 247                                           -1               ,
 248                                           Layer2TowerID(hit->GetId(),kTRUE),
 249                                           Digitize( hit->GetEnergy() ),
 250                                           hit->GetTime() ) ;
 251
 252         if((hit->GetId()/geom->GetNPhi()) < (2*geom->GetNZ()))
 253           { //cout << "ID for Preshower = " << curSDigit->GetId()  << endl ;
 254             for(Int_t check= 0; check < nSdigits; check++) {
 255               sdigit = (AliEMCALDigit *)sdigits->At(check);
 256               if( sdigit->GetId() == curSDigit->GetId())
 257                 {
 258                   *sdigit = *sdigit + *curSDigit ;
 259                   newsdigit = kFALSE ;
 260                 }
 261             }
 262             if (newsdigit)
 263               { new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
 264               nSdigits++ ;
 265               //cout << "Preshower nsdigits = " << nSdigits << endl ;
 266               }
 267             newsdigit=kTRUE;
 268           }
 269       }
 270       //    } // loop over tracks
 271       sdigits->Sort() ;
 272
 273       nSdigits = sdigits->GetEntriesFast() ;
 274       sdigits->Expand(nSdigits) ;
 275
 276       //  Int_t i ;
 277       for (i = 0 ; i < nSdigits ; i++) {
 278         AliEMCALDigit * digit = (AliEMCALDigit *) sdigits->At(i) ;
 279         digit->SetIndexInList(i) ;
 280       }
 281
 282       if(gAlice->TreeS() == 0)
 283         gAlice->MakeTree("S") ;
 284
 285
 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       TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
 300       sdigitsBranch->SetTitle(sdname);
 301       cout << " AliEMCALSDigitizer::Exec sdname " << sdname << endl ;
 302
 303       if (file) {
 304         sdigitsBranch->SetFile(file);
 305         TIter next( sdigitsBranch->GetListOfBranches());
 306         TBranch * subbr;
 307         while ((subbr=(TBranch*)next())) {
 308           subbr->SetFile(file);
 309         }
 310         cwd->cd();
 311       }
 312
 313       //second - SDigitizer
 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);
 319       if (file) {
 320         sdigitizerBranch->SetFile(file);
 321         TIter next( sdigitizerBranch->GetListOfBranches());
 322         TBranch * subbr ;
 323         while ((subbr=(TBranch*)next())) {
 324           subbr->SetFile(file);
 325         }
 326         cwd->cd();
 327         delete file;
 328       }
 329
 330       sdigitsBranch->Fill() ;
 331       sdigitizerBranch->Fill() ;
 332       gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
 333
 334       if(strstr(option,"deb"))
 335         PrintSDigits(option) ;
 336
 337   }
 338
 339   if(strstr(option,"tim")){
 340     gBenchmark->Stop("EMCALSDigitizer");
 341     cout << "AliEMCALSDigitizer:" << endl ;
 342     cout << "   took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing "
 343          <<  gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ;
 344     cout << endl ;
 345   }
 346
 347
 348 }
 349 //__________________________________________________________________
 350 void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
 351
 352   // Setting title to branch SDigits
 353
 354   TString stitle(title) ;
 355
 356   // check if branch with title already exists
 357   TBranch * sdigitsBranch    =
 358     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ;
 359   TBranch * sdigitizerBranch =
 360     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
 361   const char * sdigitsTitle    = sdigitsBranch ->GetTitle() ;
 362   const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
 363   if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
 364     cerr << "ERROR: AliEMCALSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
 365     return ;
 366   }
 367
 368   cout << "AliEMCALSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
 369
 370   SetName(title) ;
 371
 372   // Post to the WhiteBoard
 373   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 374   gime->PostSDigits( title, GetTitle()) ;
 375
 376
 377 }
 378 //__________________________________________________________________
 379 void AliEMCALSDigitizer::Print(Option_t* option)const{
 380   cout << "------------------- "<< GetName() << " -------------" << endl ;
 381   cout << "   Writing SDigitis to branch with title  " << GetName() << endl ;
 382   cout << "   with digitization parameters  A = " << fA << endl ;
 383   cout << "                                 B = " << fB << endl ;
 384   cout << "   Threshold for Primary assignment= " << fPrimThreshold << endl ;
 385   cout << "---------------------------------------------------"<<endl ;
 386
 387 }
 388 //__________________________________________________________________
 389 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const{
 390   if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
 391     return kTRUE ;
 392   else
 393     return kFALSE ;
 394 }
 395 //__________________________________________________________________
 396 void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
 397   //Prints list of digits produced at the current pass of AliEMCALDigitizer
 398
 399   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 400   TString sdname(GetName()) ;
 401   sdname.Remove(sdname.Index(GetTitle())-1) ;
 402   const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
 403
 404   cout << "AliEMCALSDigitiser: event " << gAlice->GetEvNumber() << endl ;
 405   cout << "      Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
 406   cout << endl ;
 407   if(strstr(option,"all")||strstr(option,"EMC")){
 408
 409     //loop over digits
 410     AliEMCALDigit * digit;
 411     cout << "SDigit Id " << " Amplitude " <<  " Time " << " Index "  <<  " Nprim " << " Primaries list " <<  endl;
 412     Int_t index ;
 413     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
 414       digit = (AliEMCALDigit * )  sdigits->At(index) ;
 415       cout << setw(6)  <<  digit->GetId() << "   "  <<  setw(10)  <<  digit->GetAmp() <<   "    "  << digit->GetTime()
 416            << setw(6)  <<  digit->GetIndexInList() << "    "
 417            << setw(5)  <<  digit->GetNprimary() <<"    ";
 418
 419       Int_t iprimary;
 420       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
 421         cout << " "  <<  digit->GetPrimary(iprimary+1) << "  ";
 422       cout << endl;
 423     }
 424     cout <<endl;
 425   }
 426 }
 427 //________________________________________________________________________
 428 Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t s, Bool_t preshower)
 429 {
 430   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
 431   const AliEMCALGeometry * geom = gime->EMCALGeometry();
 432   Int_t a  = (s/geom->GetNPhi())%geom->GetNZ()+1; // Phi Tower Index
 433   Int_t b = (s-1)%(geom->GetNPhi())+1; //Eta Tower Index
 434   Int_t t = -10;
 435
 436 if (a > 0 && b > 0)
 437 {
 438   t = a*b + preshower*geom->GetNPhi()*geom->GetNZ();
 439   return t;
 440 }
 441 else
 442 {cerr << " AliEMCALSDigitizer::Layer2TowerID() -- there is an error "<< endl << "Phi number = "
 443       << b << "Eta number = " << a << endl ;
 444 return t;}
 445 }