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