PHOS/AliPHOSSDigitizer.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
  17 /* $Id$ */
  18
  19 //_________________________________________________________________________
  20 // This is a TTask that makes SDigits out of Hits
  21 // The name of the TTask is also the title of the branch that will contain
  22 // the created SDigits
  23 // The title of the TTAsk is the name of the file that contains the hits from
  24 // which the SDigits are created
  25 // A Summable Digits is the sum of all hits originating
  26 // from one primary in one active cell
  27 // A threshold for assignment of the primary to SDigit is applied
  28 // SDigits are written to TreeS, branch "PHOS"
  29 // AliPHOSSDigitizer with all current parameters is written
  30 // to TreeS branch "AliPHOSSDigitizer".
  31 // Both branches have the same title. If necessary one can produce
  32 // another set of SDigits with different parameters. Two versions
  33 // can be distunguished using titles of the branches.
  34 // User case:
  35 //  root [0] AliPHOSSDigitizer * s = new AliPHOSSDigitizer("galice.root")
  36 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  37 //  root [1] s->ExecuteTask()
  38 //             // Makes SDigitis for all events stored in galice.root
  39 //  root [2] s->SetPedestalParameter(0.001)
  40 //             // One can change parameters of digitization
  41 // root [3] s->SetSDigitsBranch("Pedestal 0.001")
  42 //             // and write them into the new branch
  43 // root [4] s->ExecuteTask("deb all tim")
  44 //             // available parameters:
  45 //             deb - print # of produced SDigitis
  46 //             deb all  - print # and list of produced SDigits
  47 //             tim - print benchmarking information
  48 //
  49 //*-- Author :  Dmitri Peressounko (SUBATECH & KI)
  50 //////////////////////////////////////////////////////////////////////////////
  51
  52
  53 // --- ROOT system ---
  54 #include "TFile.h"
  55 #include "TROOT.h"
  56 #include "TBenchmark.h"
  57
  58 // --- Standard library ---
  59
  60 // --- AliRoot header files ---
  61 #include "AliRun.h"
  62 #include "AliPHOSDigit.h"
  63 #include "AliPHOSGetter.h"
  64 #include "AliPHOSHit.h"
  65 #include "AliPHOSSDigitizer.h"
  66
  67 ClassImp(AliPHOSSDigitizer)
  68
  69
  70 //____________________________________________________________________________
  71   AliPHOSSDigitizer::AliPHOSSDigitizer():TTask("","")
  72 {
  73   // ctor
  74   InitParameters() ;
  75   fDefaultInit = kTRUE ;
  76 }
  77
  78 //____________________________________________________________________________
  79 AliPHOSSDigitizer::AliPHOSSDigitizer(const char * headerFile, const char * sDigitsTitle, const Bool_t toSplit):
  80 TTask(sDigitsTitle, headerFile)
  81 {
  82   // ctor
  83   InitParameters() ;
  84   fToSplit = toSplit ;
  85   Init();
  86   fDefaultInit = kFALSE ;
  87 }
  88
  89 //____________________________________________________________________________
  90 AliPHOSSDigitizer::AliPHOSSDigitizer(const AliPHOSSDigitizer & sd) {
  91   //cpy ctor
  92
  93   fA             = sd.fA ;
  94   fB             = sd.fB ;
  95   fPrimThreshold = sd.fPrimThreshold ;
  96   fSDigitsInRun  = sd.fSDigitsInRun ;
  97   fSplitFile     = new TFile( (sd.fSplitFile)->GetName(), "new")  ;
  98   fToSplit       = sd.fToSplit ;
  99 }
 100
 101 //____________________________________________________________________________
 102 AliPHOSSDigitizer::~AliPHOSSDigitizer()
 103 {
 104   // dtor
 105
 106   fSplitFile = 0 ;
 107 }
 108
 109 //____________________________________________________________________________
 110 void AliPHOSSDigitizer::Init()
 111 {
 112   // Initialization: open root-file, allocate arrays for hits and sdigits,
 113   // attach task SDigitizer to the list of PHOS tasks
 114   //
 115   // Initialization can not be done in the default constructor
 116   //============================================================= YS
 117   //  The initialisation is now done by AliPHOSGetter
 118
 119   if( strcmp(GetTitle(), "") == 0 )
 120     SetTitle("galice.root") ;
 121
 122   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), GetName(),fToSplit) ;
 123   if ( gime == 0 ) {
 124     Error("Init" ,"Could not obtain the Getter object !") ;
 125     return ;
 126   }
 127
 128   gime->PostSDigits( GetName(), GetTitle() ) ;
 129
 130   fSplitFile = 0 ;
 131   if(fToSplit){
 132     // construct the name of the file as /path/PHOS.SDigits.root
 133     // First - extract full path if necessary
 134     TString sDigitsFileName(GetTitle()) ;
 135     Ssiz_t islash = sDigitsFileName.Last('/') ;
 136     if(islash<sDigitsFileName.Length())
 137       sDigitsFileName.Remove(islash+1,sDigitsFileName.Length()) ;
 138     else
 139       sDigitsFileName="" ;
 140     // Next - append the file name
 141     sDigitsFileName+="PHOS.SDigits." ;
 142     if((strcmp(GetName(),"Default")!=0)&&(strcmp(GetName(),"")!=0)){
 143       sDigitsFileName+=GetName() ;
 144       sDigitsFileName+="." ;
 145     }
 146     sDigitsFileName+="root" ;
 147     // Finally - check if the file already opened or open the file
 148     fSplitFile = static_cast<TFile*>(gROOT->GetFile(sDigitsFileName.Data()));
 149     if(!fSplitFile)
 150       fSplitFile =  TFile::Open(sDigitsFileName.Data(),"update") ;
 151   }
 152
 153   TString sdname(GetName() );
 154   sdname.Append(":") ;
 155   sdname.Append(GetTitle() ) ;
 156   SetName(sdname) ;
 157   gime->PostSDigitizer(this) ;
 158 }
 159
 160 //____________________________________________________________________________
 161 void AliPHOSSDigitizer::InitParameters()
 162 {
 163   // initializes the parameters for difitization
 164   fA             = 0;
 165   fB             = 10000000.;
 166   fPrimThreshold = 0.01 ;
 167   fSDigitsInRun  = 0 ;
 168   fSplitFile     = 0 ;
 169   fToSplit       = kFALSE ;
 170 }
 171
 172 //____________________________________________________________________________
 173 void AliPHOSSDigitizer::Exec(Option_t *option)
 174 {
 175   // Collects all hits in the same active volume into digit
 176
 177   if( strcmp(GetName(), "") == 0 )
 178     Init() ;
 179
 180   if (strstr(option, "print") ) {
 181     Print("") ;
 182     return ;
 183   }
 184
 185   if(strstr(option,"tim"))
 186     gBenchmark->Start("PHOSSDigitizer");
 187
 188   //Check, if this branch already exits
 189   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 190   if(gime->BranchExists("SDigits") )
 191     return;
 192
 193   TString sdname(GetName()) ;
 194   sdname.Remove(sdname.Index(GetTitle())-1) ;
 195
 196   Int_t nevents = gime->MaxEvent() ;
 197   Int_t ievent ;
 198   for(ievent = 0; ievent < nevents; ievent++){
 199     gime->Event(ievent,"H") ;
 200     const TClonesArray * hits = gime->Hits() ;
 201     TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
 202     sdigits->Clear();
 203     Int_t nSdigits = 0 ;
 204
 205     //Now make SDigits from hits, for PHOS it is the same, so just copy
 206     Int_t nPrim =  static_cast<Int_t>((gAlice->TreeH())->GetEntries()) ;
 207     // Attention nPrim is the number of primaries tracked by Geant
 208     // and this number could be different to the number of Primaries in TreeK;
 209     Int_t iprim ;
 210     for (iprim = 0 ; iprim < nPrim ; iprim ++) {
 211       //=========== Get the PHOS branch from Hits Tree for the Primary iprim
 212       gime->Track(iprim) ;
 213       Int_t i;
 214       for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
 215         AliPHOSHit * hit = dynamic_cast<AliPHOSHit *>(hits->At(i)) ;
 216         // Assign primary number only if contribution is significant
 217
 218         if( hit->GetEnergy() > fPrimThreshold)
 219           new((*sdigits)[nSdigits]) AliPHOSDigit(hit->GetPrimary(),hit->GetId(),
 220                                                  Digitize(hit->GetEnergy()), hit->GetTime()) ;
 221         else
 222           new((*sdigits)[nSdigits]) AliPHOSDigit( -1              , hit->GetId(),
 223                                                   Digitize(hit->GetEnergy()), hit->GetTime()) ;
 224         nSdigits++ ;
 225
 226       }
 227
 228     } // loop over iprim
 229
 230     sdigits->Sort() ;
 231
 232     nSdigits = sdigits->GetEntriesFast() ;
 233     fSDigitsInRun += nSdigits ;
 234     sdigits->Expand(nSdigits) ;
 235
 236     Int_t i ;
 237     for (i = 0 ; i < nSdigits ; i++) {
 238       AliPHOSDigit * digit = dynamic_cast<AliPHOSDigit *>(sdigits->At(i)) ;
 239       digit->SetIndexInList(i) ;
 240     }
 241
 242     //Now write SDigits
 243
 244     if((gAlice->TreeS() == 0)|| (fSplitFile))
 245       gAlice->MakeTree("S", fSplitFile);
 246
 247     if(fSplitFile)
 248       fSplitFile->cd() ;
 249
 250     //First list of sdigits
 251     Int_t bufferSize = 32000 ;
 252     TBranch * sdigitsBranch = gAlice->TreeS()->Branch("PHOS",&sdigits,bufferSize);
 253     sdigitsBranch->SetTitle(sdname);
 254
 255     //Next - SDigitizer
 256     Int_t splitlevel = 0 ;
 257     AliPHOSSDigitizer * sd = this ;
 258     TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliPHOSSDigitizer","AliPHOSSDigitizer",
 259                                                          &sd,bufferSize,splitlevel);
 260     sdigitizerBranch->SetTitle(sdname);
 261
 262     sdigitsBranch->Fill() ;
 263     sdigitizerBranch->Fill() ;
 264
 265     gAlice->TreeS()->AutoSave() ;
 266
 267     if(strstr(option,"deb"))
 268       PrintSDigits(option) ;
 269   }
 270
 271   if(strstr(option,"tim")){
 272     gBenchmark->Stop("PHOSSDigitizer");
 273     Info("Exec","   took %f seconds for SDigitizing  %f seconds per event",
 274          gBenchmark->GetCpuTime("PHOSSDigitizer"), gBenchmark->GetCpuTime("PHOSSDigitizer")/nevents) ;
 275   }
 276 }
 277
 278 //__________________________________________________________________
 279 void AliPHOSSDigitizer::SetSDigitsBranch(const char * title )
 280 {
 281   // Setting title to branch SDigits
 282
 283   TString stitle(title) ;
 284
 285   // check if branch with title already exists
 286   TBranch * sdigitsBranch    =
 287     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("PHOS")) ;
 288   TBranch * sdigitizerBranch =
 289     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliPHOSSDigitizer")) ;
 290   const char * sdigitsTitle    = sdigitsBranch ->GetTitle() ;
 291   const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
 292   if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
 293     Error("SetSDigitsBranch", "Cannot overwrite existing branch with title %s", title) ;
 294     return ;
 295   }
 296
 297   Info("SetSDigitsBranch", "-> Changing SDigits file from %s to %s", GetName(), title) ;
 298
 299   SetName(title) ;
 300
 301   // Post to the WhiteBoard
 302   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 303   gime->PostSDigits( title, GetTitle()) ;
 304 }
 305
 306
 307 //__________________________________________________________________
 308 void AliPHOSSDigitizer::Print(Option_t* option)const
 309 {
 310   // Prints parameters of SDigitizer
 311   TString message ;
 312   message  = "\n------------------- %s -------------\n" ;
 313   message += "   Writing SDigits to branch with title  %s\n" ;
 314   message += "   with digitization parameters  A = %f\n" ;
 315   message += "                                 B = %f\n" ;
 316   message += "   Threshold for Primary assignment= %f\n" ;
 317   message += "---------------------------------------------------\n" ;
 318   Info("Print", message.Data(),  GetName(),  GetName(), fA, fB, fPrimThreshold ) ;
 319
 320 }
 321
 322 //__________________________________________________________________
 323 Bool_t AliPHOSSDigitizer::operator==( AliPHOSSDigitizer const &sd )const
 324 {
 325   // Equal operator.
 326   // SDititizers are equal if their pedestal, slope and threshold are equal
 327
 328   if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
 329     return kTRUE ;
 330   else
 331     return kFALSE ;
 332 }
 333
 334 //__________________________________________________________________
 335 void AliPHOSSDigitizer::PrintSDigits(Option_t * option)
 336 {
 337   // Prints list of digits produced in the current pass of AliPHOSDigitizer
 338
 339
 340   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 341   TString sdname(GetName()) ;
 342   sdname.Remove(sdname.Index(GetTitle())-1) ;
 343   const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
 344
 345   TString message ;
 346   message  = "\nAliPHOSSDigitiser: event " ;
 347   message += gAlice->GetEvNumber();
 348   message += "\n      Number of entries in SDigits list " ;
 349   message += sdigits->GetEntriesFast() ;
 350   char * tempo = new char[8192];
 351
 352   if(strstr(option,"all")||strstr(option,"EMC")){
 353
 354     //loop over digits
 355     AliPHOSDigit * digit;
 356     message += "\nEMC sdigits\n" ;
 357     message += "\n   Id  Amplitude    Time          Index Nprim: Primaries list \n" ;
 358     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
 359     Int_t index ;
 360     for (index = 0 ; (index < sdigits->GetEntriesFast()) &&
 361          ((dynamic_cast<AliPHOSDigit *> (sdigits->At(index)))->GetId() <= maxEmc) ; index++) {
 362       digit = dynamic_cast<AliPHOSDigit *>( sdigits->At(index) ) ;
 363       if(digit->GetNprimary() == 0)
 364         continue;
 365       sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
 366               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 367       message += tempo ;
 368       Int_t iprimary;
 369       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 370         sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
 371         message += tempo ;
 372       }
 373     }
 374   }
 375
 376   if(strstr(option,"all")||strstr(option,"CPV")){
 377
 378     //loop over CPV digits
 379     AliPHOSDigit * digit;
 380
 381     message += "CPV sdigits\n" ;
 382     message += "\n   Id  Amplitude    Index Nprim: Primaries list \n" ;
 383     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
 384     Int_t index ;
 385     for (index = 0 ; index < sdigits->GetEntriesFast(); index++) {
 386       digit = dynamic_cast<AliPHOSDigit *>( sdigits->At(index) ) ;
 387       if(digit->GetId() > maxEmc){
 388         sprintf(tempo, "\n%6d  %8d    %4d      %2d :",
 389                 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
 390         message += tempo ;
 391         Int_t iprimary;
 392         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 393           sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
 394           message += tempo ;
 395         }
 396       }
 397     }
 398   }
 399   delete tempo ;
 400   Info("PrintSDigits", message.Data() ) ;
 401 }
 402
 403 //____________________________________________________________________________
 404 void AliPHOSSDigitizer::UseHitsFrom(const char * filename)
 405 {
 406   SetTitle(filename) ;
 407   Init() ;
 408 }