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 "TBenchmark.h"
  53 #include "TH1.h"
  54 #include "TBrowser.h"
  55 //#include "TObjectTable.h"
  56
  57 // --- Standard library ---
  58 #include "stdlib.h"
  59
  60 // --- AliRoot header files ---
  61 #include "AliLog.h"
  62 #include "AliEMCALDigit.h"
  63 #include "AliEMCALGetter.h"
  64 #include "AliEMCALHit.h"
  65 #include "AliEMCALSDigitizer.h"
  66 #include "AliEMCALGeometry.h"
  67 #include "AliEMCALJetMicroDst.h"
  68 //#include "AliMemoryWatcher.h"
  69
  70 ClassImp(AliEMCALSDigitizer)
  71
  72 //____________________________________________________________________________
  73   AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","")
  74 {
  75   // ctor
  76   Info("AliEMCALSDigitizer()", " CTOR # 1 #");
  77   fFirstEvent = fLastEvent  = fControlHists = 0;
  78   fDefaultInit = kTRUE ;
  79   fHists = 0;
  80 }
  81
  82 //____________________________________________________________________________
  83 AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName,
  84                                        const char * eventFolderName):
  85   TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
  86   fEventFolderName(eventFolderName)
  87 {
  88   // ctor
  89   Info("AliEMCALSDigitizer()", " CTOR # 2 #");
  90   fFirstEvent = fLastEvent  = fControlHists = 0 ; // runs one event by defaut
  91   Init();
  92   InitParameters() ;
  93   fDefaultInit = kFALSE ;
  94   fHists = 0;
  95 }
  96
  97
  98 //____________________________________________________________________________
  99 AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) {
 100   //cpy ctor
 101   Info("AliEMCALSDigitizer()", " CPY CTOR # 3 #");
 102
 103   fFirstEvent    = sd.fFirstEvent ;
 104   fLastEvent     = sd.fLastEvent ;
 105   fA             = sd.fA ;
 106   fB             = sd.fB ;
 107   fECPrimThreshold  = sd.fECPrimThreshold ;
 108   fSDigitsInRun  = sd.fSDigitsInRun ;
 109   SetName(sd.GetName()) ;
 110   SetTitle(sd.GetTitle()) ;
 111   fEventFolderName = sd.fEventFolderName;
 112 }
 113
 114
 115 //____________________________________________________________________________
 116 AliEMCALSDigitizer::~AliEMCALSDigitizer() {
 117   // dtor
 118   AliEMCALGetter * gime =
 119     //   AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
 120     AliEMCALGetter::Instance();
 121   gime->EmcalLoader()->CleanSDigitizer();
 122 }
 123
 124 //____________________________________________________________________________
 125 void AliEMCALSDigitizer::Init(){
 126   // Initialization: open root-file, allocate arrays for hits and sdigits,
 127   // attach task SDigitizer to the list of EMCAL tasks
 128   //
 129   // Initialization can not be done in the default constructor
 130   //============================================================= YS
 131   //  The initialisation is now done by the getter
 132
 133   fInit = kTRUE ;
 134
 135   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
 136   if ( gime == 0 ) {
 137     Fatal("Init", "Could not obtain the Getter objectfor file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;
 138     return ;
 139   }
 140
 141   TString opt("SDigits") ;
 142   if(gime->VersionExists(opt) ) {
 143     Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
 144     fInit = kFALSE ;
 145   }
 146
 147   gime->PostSDigitizer(this);
 148   gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 149
 150 }
 151
 152 //____________________________________________________________________________
 153 void AliEMCALSDigitizer::InitParameters()
 154 {
 155   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 156   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
 157   if (geom->GetSampling() == 0.) {
 158     Fatal("InitParameters", "Sampling factor not set !") ;
 159   }
 160   // Initializes parameters
 161   fA         = 0;
 162   fB         = 1.e+7;
 163   fSampling  = geom->GetSampling();
 164
 165  // threshold for deposit energy of hit
 166   fECPrimThreshold  = 0.; // 24-nov-04 - was 1.e-6;
 167 }
 168
 169 //____________________________________________________________________________
 170 void AliEMCALSDigitizer::Exec(Option_t *option)
 171 {
 172   // Collects all hit of the same tower into digits
 173   TString o(option); o.ToUpper();
 174   if (strstr(option, "print") ) {
 175     Print() ;
 176     return ;
 177   }
 178
 179   if(strstr(option,"tim"))
 180     gBenchmark->Start("EMCALSDigitizer");
 181
 182   //AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ;
 183   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 184
 185   //switch off reloading of this task while getting event
 186   if (!fInit) { // to prevent overwrite existing file
 187     AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
 188     return ;
 189     }
 190
 191   if (fLastEvent == -1)
 192     fLastEvent = gime->MaxEvent() - 1 ;
 193   else {
 194     if(o != "EXACT") fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent());
 195     fLastEvent = TMath::Min(fLastEvent, gime->MaxEvent());
 196     printf("AliEMCALSDigitizer::Exec : option: %s | %i -> %i events : Max events %i \n \n",
 197     option, fFirstEvent, fLastEvent, gime->MaxEvent());
 198   }
 199   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 200
 201   Int_t ievent;
 202   Float_t energy=0.; // de * fSampling - 23-nov-04
 203   for (ievent = fFirstEvent; ievent < fLastEvent; ievent++) {
 204     if(ievent%100==0 || ievent==fLastEvent-1) printf(" processed event %i \n", ievent);
 205     gime->Event(ievent,"XH"); // read primaries and hits onl
 206     TTree * treeS = gime->TreeS();
 207     TClonesArray * hits = gime->Hits() ;
 208     TClonesArray * sdigits = gime->SDigits() ;
 209     sdigits->Clear();
 210     Int_t nSdigits = 0 ;
 211     //Now make SDigits from hits, for EMCAL it is the same, so just copy
 212     Int_t nPrim =  static_cast<Int_t>((gime->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;
 215     Int_t iprim ;
 216     for ( iprim = 0 ; iprim < nPrim ; iprim++ ) {
 217       //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
 218       gime->Track(iprim) ;
 219       Int_t i;
 220       AliEMCALGeometry *geom = gime->EMCALGeometry();
 221       for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
 222         AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
 223         AliEMCALDigit * curSDigit = 0 ;
 224         AliEMCALDigit * sdigit = 0 ;
 225         Bool_t newsdigit = kTRUE;
 226
 227         // hit->GetId() - Absolute Id number EMCAL segment
 228         if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId())
 229           energy = hit->GetEnergy() * fSampling; // 23-nov-04
 230           if(energy >  fECPrimThreshold )
 231         // Assign primary number only if deposited energy is significant
 232             curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
 233                                             hit->GetIparent(), hit->GetId(),
 234                                             Digitize(energy), hit->GetTime() ) ;
 235           else
 236             curSDigit =  new AliEMCALDigit( -1               ,
 237                                             -1               ,
 238                                             hit->GetId(),
 239                                             Digitize(energy), hit->GetTime() ) ;
 240         } else {
 241           Warning("Exec"," abs id %i is bad \n", hit->GetId());
 242           newsdigit = kFALSE;
 243         }
 244
 245         Int_t check = 0 ;
 246
 247         for(check= 0; check < nSdigits ; check++) {
 248           sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
 249           if(curSDigit != 0){
 250             if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?
 251             *sdigit = *sdigit + *curSDigit;
 252             newsdigit = kFALSE;
 253             }
 254           }
 255         }
 256         if (newsdigit) {
 257           new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
 258           nSdigits++ ;
 259         }
 260         delete curSDigit ;
 261       }  // loop over all hits (hit = deposited energy/entering particle)
 262     } // loop over iprim
 263     sdigits->Sort() ;
 264
 265     nSdigits = sdigits->GetEntriesFast() ;
 266     fSDigitsInRun += nSdigits ;
 267     sdigits->Expand(nSdigits) ;
 268
 269     Int_t nPrimarymax = -1 ;
 270     Int_t i ;
 271     Double_t e=0.,esum=0.;
 272     sv::FillH1(fHists, 0, double(sdigits->GetEntriesFast()));
 273     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {
 274       AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
 275       sdigit->SetIndexInList(i) ;
 276
 277       sv::FillH1(fHists, 2, double(sdigit->GetAmp()));
 278       e = double(Calibrate(sdigit->GetAmp()));
 279       esum += e;
 280       sv::FillH1(fHists, 3, e);
 281       sv::FillH1(fHists, 4, double(sdigit->GetId()));
 282     }
 283     if(esum>0.) sv::FillH1(fHists, 1, esum);
 284
 285     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { // for what
 286       if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > nPrimarymax)
 287         nPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
 288     }
 289
 290     // Now write SDigits
 291     //First list of sdigits
 292
 293     Int_t bufferSize = 32000 ;
 294     TBranch * sdigitsBranch = treeS->Branch("EMCAL",&sdigits,bufferSize);
 295
 296     sdigitsBranch->Fill() ;
 297
 298     gime->WriteSDigits("OVERWRITE");
 299
 300     //NEXT - SDigitizer
 301     if(ievent == fFirstEvent) gime->WriteSDigitizer("OVERWRITE");  // why in event cycle ?
 302
 303     if(strstr(option,"deb"))
 304       PrintSDigits(option) ;
 305   }
 306   //  gime->WriteSDigitizer("OVERWRITE"); // 12-jan-04
 307
 308   Unload();
 309
 310   gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
 311
 312   if(strstr(option,"tim")){
 313     gBenchmark->Stop("EMCALSDigitizer");
 314     printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n",
 315          gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ;
 316   }
 317
 318   //TFile f("out.root","RECREATE");
 319   //memwatcher.WriteToFile();
 320   //f.Close();
 321 }
 322
 323
 324 void AliEMCALSDigitizer::Print1(Option_t * option)
 325 {
 326   Print();
 327   PrintSDigits(option);
 328 }
 329
 330 //__________________________________________________________________
 331 void AliEMCALSDigitizer::Print() const
 332 {
 333   // Prints parameters of SDigitizer
 334   printf("Print: \n------------------- %s -------------\n", GetName() ) ;
 335   printf("   fInit                                 %i\n", int(fInit));
 336   printf("   fFirstEvent                           %i\n", fFirstEvent);
 337   printf("   fLastEvent                            %i\n", fLastEvent);
 338   printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
 339   printf("   with digitization parameters       A = %f\n", fA) ;
 340   printf("                                      B = %f\n", fB) ;
 341   printf("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold)  ;
 342   printf("   Sampling                             = %f\n", fSampling);
 343   printf("---------------------------------------------------\n") ;
 344 }
 345
 346 //__________________________________________________________________
 347 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
 348 {
 349   // Equal operator.
 350   // SDititizers are equal if their pedestal, slope and threshold are equal
 351   if( (fA==sd.fA)&&(fB==sd.fB)&&
 352       (fECPrimThreshold==sd.fECPrimThreshold))
 353     return kTRUE ;
 354   else
 355     return kFALSE ;
 356 }
 357
 358 //__________________________________________________________________
 359 void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
 360 {
 361   //Prints list of digits produced at the current pass of AliEMCALDigitizer
 362
 363
 364   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 365   const TClonesArray * sdigits = gime->SDigits() ;
 366
 367   printf("\n") ;
 368   printf("event %i", gAlice->GetEvNumber()) ;
 369   printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast());
 370   if(strstr(option,"all")||strstr(option,"EMC")){
 371
 372     //loop over digits
 373     AliEMCALDigit * digit;
 374     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;
 375     Int_t index, isum=0;
 376     char * tempo = new char[8192];
 377     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
 378       digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
 379       sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
 380               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
 381       printf(tempo);
 382       isum += digit->GetAmp();
 383
 384       Int_t iprimary;
 385       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
 386         sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ;
 387         printf(tempo);
 388       }
 389     }
 390     delete tempo ;
 391     printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6);
 392   } else printf("\n");
 393 }
 394
 395 //____________________________________________________________________________
 396 void AliEMCALSDigitizer::Unload() const
 397 {
 398   // Unload Hits and SDigits from the folder
 399   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 400   AliEMCALLoader * loader = gime->EmcalLoader() ;
 401   loader->UnloadHits() ;
 402   loader->UnloadSDigits() ;
 403 }
 404
 405 void AliEMCALSDigitizer::Browse(TBrowser* b)
 406 {
 407   if(fHists) b->Add(fHists);
 408   TTask::Browse(b);
 409 }
 410
 411 TList *AliEMCALSDigitizer::BookControlHists(const int var)
 412 { // 22-nov-04
 413   gROOT->cd();
 414   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName);
 415   const AliEMCALGeometry *geom = gime->EMCALGeometry() ;
 416   if(var>=1){
 417     new TH1F("HSDigiN",  "#EMCAL  sdigits ", 1001, -0.5, 1000.5);
 418     new TH1F("HSDigiSumEnergy","Sum.EMCAL energy", 1000, 0.0, 100.);
 419     new TH1F("HSDigiAmp",  "EMCAL sdigits amplitude", 1000, 0., 2.e+9);
 420     new TH1F("HSDigiEnergy","EMCAL cell energy", 1000, 0.0, 100.);
 421     new TH1F("HSDigiAbsId","EMCAL absID for sdigits",
 422     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
 423   }
 424   fHists = sv::MoveHistsToList("EmcalSDigiControlHists", kFALSE);
 425   return fHists;
 426 }
 427
 428 void AliEMCALSDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
 429 {
 430   sv::SaveListOfHists(fHists, name, kSingleKey, opt);
 431 }