PHOS/AliPHOSGetter.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 //  A singleton. This class should be used in the analysis stage to get
  20 //  reconstructed objects: Digits, RecPoints, TrackSegments and RecParticles,
  21 //  instead of directly reading them from galice.root file. This container
  22 //  ensures, that one reads Digits, made of these particular digits, RecPoints,
  23 //  made of these particular RecPoints, TrackSegments and RecParticles.
  24 //  This becomes non trivial if there are several identical branches, produced with
  25 //  different set of parameters.
  26 //
  27 //  An example of how to use (see also class AliPHOSAnalyser):
  28 //  AliPHOSGetter * gime = AliPHOSGetter::GetInstance("galice.root","test") ;
  29 //  for(Int_t irecp = 0; irecp < gime->NRecParticles() ; irecp++)
  30 //     AliPHOSRecParticle * part = gime->RecParticle(1) ;
  31 //     ................
  32 //  gime->Event(event) ;    // reads new event from galice.root
  33 //
  34 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC KI & SUBATECH)
  35 //*--         Completely redesigned by Dmitri Peressounko March 2001
  36 //
  37 //*-- YS June 2001 : renamed the original AliPHOSIndexToObject and make
  38 //*--         systematic usage of TFolders without changing the interface
  39 //////////////////////////////////////////////////////////////////////////////
  40
  41 // --- ROOT system ---
  42
  43 #include <TFile.h>
  44 #include <TROOT.h>
  45 #include <TSystem.h>
  46 #include <TParticle.h>
  47 #include <TF1.h>
  48 #include <TGraph.h>
  49 #include <TCanvas.h>
  50 #include <TStyle.h>
  51 //#include <TFrame.h>
  52
  53 // --- Standard library ---
  54
  55 // --- AliRoot header files ---
  56 #include "AliESD.h"
  57 #include "AliHeader.h"
  58 #include "AliMC.h"
  59 #include "AliPHOS.h"
  60 #include "AliPHOSBeamTestEvent.h"
  61 #include "AliPHOSGetter.h"
  62 #include "AliPHOSLoader.h"
  63 #include "AliPHOSPulseGenerator.h"
  64 #include "AliRunLoader.h"
  65 #include "AliStack.h"
  66 #include "AliCaloRawStream.h"
  67 #include "AliRawReaderFile.h"
  68 #include "AliLog.h"
  69 #include "AliCDBLocal.h"
  70 #include "AliCDBStorage.h"
  71 #include "AliCDBManager.h"
  72
  73 ClassImp(AliPHOSGetter)
  74
  75 AliPHOSGetter   * AliPHOSGetter::fgObjGetter  = 0 ;
  76 AliPHOSLoader   * AliPHOSGetter::fgPhosLoader = 0;
  77 AliPHOSCalibData* AliPHOSGetter::fgCalibData  = 0;
  78 Int_t AliPHOSGetter::fgDebug = 0;
  79
  80 //  TFile * AliPHOSGetter::fgFile = 0 ;
  81
  82 AliPHOSGetter::AliPHOSGetter() :
  83   fBTE(0),
  84   fLoadingStatus(),
  85   fNPrimaries(0),
  86   fPrimaries(0),
  87   fESDFile(0),
  88   fESDFileName(),
  89   fESD(0),
  90   fESDTree(0),
  91   fRawDigits(kFALSE),
  92   fcdb(0)
  93 {
  94   // ctor: this is a singleton, the ctor should never be called but cint needs it as public
  95   Fatal("ctor", "AliPHOSGetter is a singleton default ctor not callable") ;
  96 }
  97
  98
  99 //____________________________________________________________________________
 100 AliPHOSGetter::AliPHOSGetter(const char* headerFile, const char* version, Option_t * openingOption) :
 101   fBTE(0),
 102   fLoadingStatus(),
 103   fNPrimaries(0),
 104   fPrimaries(0),
 105   fESDFile(0),
 106   fESDFileName(),
 107   fESD(0),
 108   fESDTree(0),
 109   fRawDigits(kFALSE),
 110   fcdb(0)
 111 {
 112   // ctor only called by Instance()
 113
 114   AliRunLoader* rl = AliRunLoader::GetRunLoader(version) ;
 115   if (!rl) {
 116     rl = AliRunLoader::Open(headerFile, version, openingOption);
 117     if (!rl) {
 118       Fatal("AliPHOSGetter", "Could not find the Run Loader for %s - %s",headerFile, version) ;
 119       return ;
 120     }
 121     if (rl->GetAliRun() == 0x0) {
 122       rl->LoadgAlice();
 123       gAlice = rl->GetAliRun(); // should be removed
 124       rl->LoadHeader();
 125     }
 126   }
 127   fgPhosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 128   if ( !fgPhosLoader )
 129     Error("AliPHOSGetter", "Could not find PHOSLoader") ;
 130   else
 131     fgPhosLoader->SetTitle(version);
 132
 133   // initialize data members
 134   SetDebug(0) ;
 135   fBTE = 0 ;
 136   fPrimaries = 0 ;
 137   fLoadingStatus = "" ;
 138
 139   fESDFileName = rl->GetFileName()  ; // this should be the galice.root file
 140   fESDFileName.ReplaceAll("galice.root", "AliESDs.root") ;
 141   fESDFile = 0 ;
 142   fESD = 0 ;
 143   fESDTree = 0 ;
 144   fRawDigits = kFALSE ;
 145
 146 }
 147
 148 AliPHOSGetter::AliPHOSGetter(const AliPHOSGetter & obj) :
 149   TObject(obj),
 150   fBTE(0),
 151   fLoadingStatus(),
 152   fNPrimaries(0),
 153   fPrimaries(0),
 154   fESDFile(0),
 155   fESDFileName(),
 156   fESD(0),
 157   fESDTree(0),
 158   fRawDigits(kFALSE),
 159   fcdb(0)
 160 {
 161   // cpy ctor requested by Coding Convention
 162   Fatal("cpy ctor", "not implemented") ;
 163 }
 164
 165 //____________________________________________________________________________
 166 AliPHOSGetter::AliPHOSGetter(Int_t /*i*/) :
 167   fBTE(0),
 168   fLoadingStatus(),
 169   fNPrimaries(0),
 170   fPrimaries(0),
 171   fESDFile(0),
 172   fESDFileName(),
 173   fESD(0),
 174   fESDTree(0),
 175   fRawDigits(kFALSE),
 176   fcdb(0)
 177 {
 178   // special constructor for onflight
 179 }
 180
 181
 182 //____________________________________________________________________________
 183 AliPHOSGetter::~AliPHOSGetter()
 184 {
 185   // dtor
 186   if(fgPhosLoader){
 187     delete fgPhosLoader ;
 188     fgPhosLoader = 0 ;
 189   }
 190   if(fBTE){
 191     delete fBTE ;
 192     fBTE = 0 ;
 193   }
 194   if(fPrimaries){
 195     fPrimaries->Delete() ;
 196     delete fPrimaries ;
 197   }
 198   if (fESD)
 199     delete fESD ;
 200   if (fESDTree)
 201     delete fESDTree ;
 202
 203   fgObjGetter = 0;
 204 }
 205
 206 //____________________________________________________________________________
 207 void AliPHOSGetter::Reset()
 208 {
 209   // resets things in case the getter is called consecutively with different files
 210   // the PHOS Loader is already deleted by the Run Loader
 211
 212   if (fPrimaries) {
 213     fPrimaries->Delete() ;
 214     delete fPrimaries ;
 215   }
 216   fgPhosLoader = 0;
 217   fgObjGetter = 0;
 218 }
 219
 220 //____________________________________________________________________________
 221 AliPHOSClusterizer * AliPHOSGetter::Clusterizer()
 222 {
 223   // Returns pointer to the Clusterizer task
 224   AliPHOSClusterizer * rv ;
 225   rv =  dynamic_cast<AliPHOSClusterizer *>(PhosLoader()->Reconstructioner()) ;
 226   if (!rv) {
 227     Event(0, "R") ;
 228     rv =  dynamic_cast<AliPHOSClusterizer*>(PhosLoader()->Reconstructioner()) ;
 229   }
 230   return rv ;
 231 }
 232
 233 //____________________________________________________________________________
 234 TObjArray * AliPHOSGetter::CpvRecPoints() const
 235 {
 236   // asks the Loader to return the CPV RecPoints container
 237
 238   TObjArray * rv = 0 ;
 239
 240   rv = PhosLoader()->CpvRecPoints() ;
 241   if (!rv) {
 242     PhosLoader()->MakeRecPointsArray() ;
 243     rv = PhosLoader()->CpvRecPoints() ;
 244   }
 245   return rv ;
 246 }
 247
 248 //____________________________________________________________________________
 249 TClonesArray * AliPHOSGetter::Digits() const
 250 {
 251   // asks the Loader to return the Digits container
 252
 253   TClonesArray * rv = 0 ;
 254   rv = PhosLoader()->Digits() ;
 255
 256   if( !rv ) {
 257     PhosLoader()->MakeDigitsArray() ;
 258     rv = PhosLoader()->Digits() ;
 259   }
 260   return rv ;
 261 }
 262
 263 //____________________________________________________________________________
 264 AliPHOSDigitizer * AliPHOSGetter::Digitizer()
 265 {
 266   // Returns pointer to the Digitizer task
 267   AliPHOSDigitizer * rv ;
 268   rv =  dynamic_cast<AliPHOSDigitizer *>(PhosLoader()->Digitizer()) ;
 269   if (!rv) {
 270     Event(0, "D") ;
 271     rv =  dynamic_cast<AliPHOSDigitizer *>(PhosLoader()->Digitizer()) ;
 272   }
 273   return rv ;
 274 }
 275
 276
 277 //____________________________________________________________________________
 278 TObjArray * AliPHOSGetter::EmcRecPoints() const
 279 {
 280   // asks the Loader to return the EMC RecPoints container
 281
 282   TObjArray * rv = 0 ;
 283
 284   rv = PhosLoader()->EmcRecPoints() ;
 285   if (!rv) {
 286     PhosLoader()->MakeRecPointsArray() ;
 287     rv = PhosLoader()->EmcRecPoints() ;
 288   }
 289   return rv ;
 290 }
 291
 292 //____________________________________________________________________________
 293 TClonesArray * AliPHOSGetter::TrackSegments() const
 294 {
 295   // asks the Loader to return the TrackSegments container
 296
 297   TClonesArray * rv = 0 ;
 298
 299   rv = PhosLoader()->TrackSegments() ;
 300   if (!rv) {
 301     PhosLoader()->MakeTrackSegmentsArray() ;
 302     rv = PhosLoader()->TrackSegments() ;
 303   }
 304   return rv ;
 305 }
 306
 307 //____________________________________________________________________________
 308 AliPHOSTrackSegmentMaker * AliPHOSGetter::TrackSegmentMaker()
 309 {
 310   // Returns pointer to the TrackSegmentMaker task
 311   AliPHOSTrackSegmentMaker * rv ;
 312   rv =  dynamic_cast<AliPHOSTrackSegmentMaker *>(PhosLoader()->TrackSegmentMaker()) ;
 313   if (!rv) {
 314     Event(0, "T") ;
 315     rv =  dynamic_cast<AliPHOSTrackSegmentMaker *>(PhosLoader()->TrackSegmentMaker()) ;
 316   }
 317   return rv ;
 318 }
 319
 320 //____________________________________________________________________________
 321 TClonesArray * AliPHOSGetter::RecParticles() const
 322 {
 323   // asks the Loader to return the TrackSegments container
 324
 325   TClonesArray * rv = 0 ;
 326
 327   rv = PhosLoader()->RecParticles() ;
 328   if (!rv) {
 329     PhosLoader()->MakeRecParticlesArray() ;
 330     rv = PhosLoader()->RecParticles() ;
 331   }
 332   return rv ;
 333 }
 334 //____________________________________________________________________________
 335 void AliPHOSGetter::Event(Int_t event, const char* opt)
 336 {
 337   // Reads the content of all Tree's S, D and R
 338
 339 //   if ( event >= MaxEvent() ) {
 340 //     Error("Event", "%d not found in TreeE !", event) ;
 341 //     return ;
 342 //   }
 343
 344   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 345
 346 //   // checks if we are dealing with test-beam data
 347 //   TBranch * btb = rl->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
 348 //   if(btb){
 349 //     if(!fBTE)
 350 //       fBTE = new AliPHOSBeamTestEvent() ;
 351 //     btb->SetAddress(&fBTE) ;
 352 //     btb->GetEntry(event) ;
 353 //   }
 354 //   else{
 355 //     if(fBTE){
 356 //       delete fBTE ;
 357 //       fBTE = 0 ;
 358 //     }
 359 //   }
 360
 361   // Loads the type of object(s) requested
 362
 363   rl->GetEvent(event) ;
 364
 365   if(strstr(opt,"X") || (strcmp(opt,"")==0)){
 366     ReadPrimaries() ;
 367   }
 368
 369   if(strstr(opt,"H")  ){
 370     ReadTreeH();
 371   }
 372
 373   if(strstr(opt,"S")  ){
 374     ReadTreeS() ;
 375   }
 376
 377   if(strstr(opt,"D") ){
 378     ReadTreeD() ;
 379   }
 380
 381   if(strstr(opt,"R") ){
 382     ReadTreeR() ;
 383   }
 384
 385   if( strstr(opt,"T") ){
 386     ReadTreeT() ;
 387   }
 388
 389   if( strstr(opt,"P") ){
 390     ReadTreeP() ;
 391   }
 392
 393   if( strstr(opt,"E") ){
 394     ReadTreeE(event) ;
 395   }
 396
 397 }
 398
 399
 400 //____________________________________________________________________________
 401 void AliPHOSGetter::Event(AliRawReader *rawReader, const char* opt, Bool_t isOldRCUFormat)
 402 {
 403   // Reads the raw event from rawReader
 404   // isOldRCUFormat defines whenever to assume
 405   // the old RCU format or not
 406
 407   if( strstr(opt,"W")  ){
 408     rawReader->NextEvent();
 409     ReadRaw(rawReader,isOldRCUFormat) ;
 410   }
 411
 412 }
 413
 414
 415 //____________________________________________________________________________
 416 Int_t AliPHOSGetter::EventNumber() const
 417   {
 418   // return the current event number
 419   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 420   return static_cast<Int_t>(rl->GetEventNumber()) ;
 421 }
 422
 423 //____________________________________________________________________________
 424   TClonesArray * AliPHOSGetter::Hits() const
 425 {
 426   // asks the loader to return  the Hits container
 427
 428   TClonesArray * rv = 0 ;
 429
 430   rv = PhosLoader()->Hits() ;
 431   if ( !rv ) {
 432     PhosLoader()->LoadHits("read");
 433     rv = PhosLoader()->Hits() ;
 434   }
 435   return rv ;
 436 }
 437
 438 //____________________________________________________________________________
 439 AliPHOSGetter * AliPHOSGetter::Instance(const char* alirunFileName, const char* version, Option_t * openingOption)
 440 {
 441   // Creates and returns the pointer of the unique instance
 442   // Must be called only when the environment has changed
 443
 444   if(!fgObjGetter){ // first time the getter is called
 445     fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 446   }
 447   else { // the getter has been called previously
 448     AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle());
 449     if ( rl->GetFileName() == alirunFileName ) {// the alirunFile has the same name
 450       // check if the file is already open
 451       TFile * galiceFile = dynamic_cast<TFile *>(gROOT->FindObject(rl->GetFileName()) ) ;
 452
 453       if ( !galiceFile )
 454         fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 455
 456       else {  // the file is already open check the version name
 457         TString currentVersionName = rl->GetEventFolder()->GetName() ;
 458         TString newVersionName(version) ;
 459         if (currentVersionName == newVersionName)
 460           if(fgDebug)
 461             ::Warning( "Instance", "Files with version %s already open", currentVersionName.Data() ) ;
 462         else {
 463           fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 464         }
 465       }
 466     }
 467     else {
 468       AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle()) ;
 469       if ( strstr(version, AliConfig::GetDefaultEventFolderName()) ) // false in case of merging
 470         delete rl ;
 471       fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 472     }
 473   }
 474   if (!fgObjGetter)
 475     ::Error("AliPHOSGetter::Instance", "Failed to create the PHOS Getter object") ;
 476   else
 477     if (fgDebug)
 478       Print() ;
 479
 480   return fgObjGetter ;
 481 }
 482
 483 //____________________________________________________________________________
 484 AliPHOSGetter *  AliPHOSGetter::Instance()
 485 {
 486   // Returns the pointer of the unique instance already defined
 487
 488   if(!fgObjGetter && fgDebug)
 489      ::Warning("AliPHOSGetter::Instance", "Getter not initialized") ;
 490
 491    return fgObjGetter ;
 492
 493 }
 494
 495 //____________________________________________________________________________
 496 Int_t AliPHOSGetter::MaxEvent() const
 497 {
 498   // returns the number of events in the run (from TE)
 499
 500   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 501   return static_cast<Int_t>(rl->GetNumberOfEvents()) ;
 502 }
 503
 504 //____________________________________________________________________________
 505 TParticle * AliPHOSGetter::Primary(Int_t index) const
 506 {
 507   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 508   return rl->Stack()->Particle(index) ;
 509 }
 510
 511 //____________________________________________________________________________
 512 AliPHOS * AliPHOSGetter:: PHOS() const
 513 {
 514   // returns the PHOS object
 515   AliPHOSLoader *    loader = 0;
 516   static AliPHOSLoader * oldloader = 0;
 517   static AliPHOS * phos = 0;
 518
 519   loader = PhosLoader();
 520
 521   if ( loader != oldloader) {
 522     phos = dynamic_cast<AliPHOS*>(loader->GetModulesFolder()->FindObject("PHOS")) ;
 523     oldloader = loader;
 524   }
 525   if (!phos)
 526     if (fgDebug)
 527       Warning("PHOS", "PHOS module not found in module folders: %s", PhosLoader()->GetModulesFolder()->GetName() ) ;
 528   return phos ;
 529 }
 530
 531
 532
 533 //____________________________________________________________________________
 534 AliPHOSPID * AliPHOSGetter::PID()
 535 {
 536   // Returns pointer to the PID task
 537   AliPHOSPID * rv ;
 538   rv =  dynamic_cast<AliPHOSPID *>(PhosLoader()->PIDTask()) ;
 539   if (!rv) {
 540     Event(0, "P") ;
 541     rv =  dynamic_cast<AliPHOSPID *>(PhosLoader()->PIDTask()) ;
 542   }
 543   return rv ;
 544 }
 545
 546 //____________________________________________________________________________
 547 AliPHOSGeometry * AliPHOSGetter::PHOSGeometry() const
 548 {
 549   // Returns PHOS geometry
 550
 551   AliPHOSGeometry * rv = 0 ;
 552   if (PHOS() )
 553     rv =  PHOS()->GetGeometry() ;
 554   else {
 555     rv = AliPHOSGeometry::GetInstance();
 556     if (!rv) {
 557       AliError("Could not find PHOS geometry! Loading the default one !");
 558       rv = AliPHOSGeometry::GetInstance("IHEP","");
 559     }
 560   }
 561   return rv ;
 562 }
 563
 564 //____________________________________________________________________________
 565 TClonesArray * AliPHOSGetter::Primaries()
 566 {
 567   // creates the Primaries container if needed
 568   if ( !fPrimaries ) {
 569     if (fgDebug)
 570       Info("Primaries", "Creating a new TClonesArray for primaries") ;
 571     fPrimaries = new TClonesArray("TParticle", 1000) ;
 572   }
 573   return fPrimaries ;
 574 }
 575
 576 //____________________________________________________________________________
 577 void  AliPHOSGetter::Print()
 578 {
 579   // Print usefull information about the getter
 580
 581   AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle());
 582   ::Info( "Print", "gAlice file is %s -- version name is %s", (rl->GetFileName()).Data(), rl->GetEventFolder()->GetName() ) ;
 583 }
 584
 585 //____________________________________________________________________________
 586 void AliPHOSGetter::ReadPrimaries()
 587 {
 588   // Read Primaries from Kinematics.root
 589
 590   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 591
 592   // gets kine tree from the root file (Kinematics.root)
 593   if ( ! rl->TreeK() ) { // load treeK the first time
 594     rl->LoadKinematics() ;
 595   }
 596
 597   fNPrimaries = (rl->GetHeader())->GetNtrack();
 598   if (fgDebug)
 599     Info( "ReadTreeK", "Found %d particles in event # %d", fNPrimaries, EventNumber() ) ;
 600
 601
 602   // first time creates the container
 603   if ( Primaries() )
 604     fPrimaries->Clear() ;
 605
 606   Int_t index = 0 ;
 607   for (index = 0 ; index < fNPrimaries; index++) {
 608     new ((*fPrimaries)[index]) TParticle(*(Primary(index)));
 609   }
 610 }
 611
 612 //____________________________________________________________________________
 613 Bool_t AliPHOSGetter::OpenESDFile()
 614 {
 615   //Open the ESD file
 616   Bool_t rv = kTRUE ;
 617   if (!fESDFile) {
 618     fESDFile = TFile::Open(fESDFileName) ;
 619     if (!fESDFile )
 620       return kFALSE ;
 621   }
 622   else if (fESDFile->IsOpen()) {
 623     fESDFile->Close() ;
 624     fESDFile = TFile::Open(fESDFileName) ;
 625   }
 626   if (!fESDFile->IsOpen())
 627     rv = kFALSE ;
 628   return rv ;
 629 }
 630
 631 //____________________________________________________________________________
 632 void AliPHOSGetter::FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Double_t & energy, Double_t & time) const
 633 {
 634   // Fits the raw signal time distribution
 635
 636   const Int_t kNoiseThreshold = 0 ;
 637   Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
 638   Double_t signal = 0., signalmax = 0. ;
 639   time   = 0. ;
 640   energy = 0. ;
 641
 642   // Create a shaper pulse object which contains all the shaper parameters
 643   AliPHOSPulseGenerator pulse;
 644
 645   if (lowGainFlag) {
 646     timezero1 = timezero2 = signalmax = timemax = 0. ;
 647     signalF->FixParameter(0, pulse.GetRawFormatLowCharge()) ;
 648     signalF->FixParameter(1, pulse.GetRawFormatLowGain()) ;
 649     Int_t index ;
 650     for (index = 0; index < pulse.GetRawFormatTimeBins(); index++) {
 651       gLowGain->GetPoint(index, time, signal) ;
 652       if (signal > kNoiseThreshold && timezero1 == 0.)
 653         timezero1 = time ;
 654       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
 655         timezero2 = time ;
 656       if (signal > signalmax) {
 657         signalmax = signal ;
 658         timemax   = time ;
 659       }
 660     }
 661     signalmax /=
 662       pulse.RawResponseFunctionMax(pulse.GetRawFormatLowCharge(),
 663                                    pulse.GetRawFormatLowGain()) ;
 664     if ( timezero1 + pulse.GetRawFormatTimePeak() < pulse.GetRawFormatTimeMax() * 0.4 ) { // else its noise
 665       signalF->SetParameter(2, signalmax) ;
 666       signalF->SetParameter(3, timezero1) ;
 667       gLowGain->Fit(signalF, "QRO", "", 0., timezero2); //, "QRON") ;
 668       energy = signalF->GetParameter(2) ;
 669       time   = signalF->GetMaximumX() - pulse.GetRawFormatTimePeak() - pulse.GetRawFormatTimeTrigger() ;
 670     }
 671   } else {
 672     timezero1 = timezero2 = signalmax = timemax = 0. ;
 673     signalF->FixParameter(0, pulse.GetRawFormatHighCharge()) ;
 674     signalF->FixParameter(1, pulse.GetRawFormatHighGain()) ;
 675     Int_t index ;
 676     for (index = 0; index < pulse.GetRawFormatTimeBins(); index++) {
 677       gHighGain->GetPoint(index, time, signal) ;
 678       if (signal > kNoiseThreshold && timezero1 == 0.)
 679         timezero1 = time ;
 680       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
 681         timezero2 = time ;
 682       if (signal > signalmax) {
 683         signalmax = signal ;
 684         timemax   = time ;
 685       }
 686     }
 687     signalmax /= pulse.RawResponseFunctionMax(pulse.GetRawFormatHighCharge(),
 688                                               pulse.GetRawFormatHighGain()) ;;
 689     if ( timezero1 + pulse.GetRawFormatTimePeak() < pulse.GetRawFormatTimeMax() * 0.4 ) { // else its noise
 690       signalF->SetParameter(2, signalmax) ;
 691       signalF->SetParameter(3, timezero1) ;
 692       gHighGain->Fit(signalF, "QRO", "", 0., timezero2) ;
 693       energy = signalF->GetParameter(2) ;
 694       time   = signalF->GetMaximumX() - pulse.GetRawFormatTimePeak() - pulse.GetRawFormatTimeTrigger() ;
 695     }
 696   }
 697   if (time == 0) energy = 0 ;
 698 }
 699
 700 //____________________________________________________________________________
 701 Int_t AliPHOSGetter::CalibrateRaw(Double_t energy, Int_t *relId)
 702 {
 703   // Convert energy into digitized amplitude for a cell relId
 704   // It is a user responsilibity to open CDB and set
 705   // AliPHOSCalibData object by the following operators:
 706   //
 707   // AliCDBLocal *loc = new AliCDBLocal("deCalibDB");
 708   // AliPHOSCalibData* clb = (AliPHOSCalibData*)AliCDBStorage::Instance()
 709   //    ->Get(path_to_calibdata,run_number);
 710   // AliPHOSGetter* gime = AliPHOSGetter::Instance("galice.root");
 711   // gime->SetCalibData(clb);
 712
 713   if (CalibData() == 0)
 714     Warning("CalibrateRaw","Calibration DB is not initiated!");
 715
 716   Int_t   module = relId[0];
 717   Int_t   column = relId[3];
 718   Int_t   row    = relId[2];
 719
 720   Float_t gainFactor = 0.0015; // width of one Emc ADC channel in GeV
 721   Float_t pedestal   = 0.005;  // Emc pedestals
 722
 723   if(CalibData()) {
 724     gainFactor = CalibData()->GetADCchannelEmc (module,column,row);
 725     pedestal   = CalibData()->GetADCpedestalEmc(module,column,row);
 726   }
 727
 728   Int_t   amp = static_cast<Int_t>( (energy - pedestal) / gainFactor + 0.5 ) ;
 729   return amp;
 730 }
 731 //____________________________________________________________________________
 732 Int_t AliPHOSGetter::ReadRaw(AliRawReader *rawReader,Bool_t isOldRCUFormat)
 733 {
 734   // reads the raw format data, converts it into digits format and store digits in Digits()
 735   // container.
 736   // isOldRCUFormat = kTRUE in case of the old RCU
 737   // format used in the raw data readout.
 738   // Reimplemented by Boris Polichtchouk (Jul 2006)
 739   // to make it working with the Jul-Aug 2006 beam test data.
 740
 741   AliCaloRawStream in(rawReader,"PHOS");
 742   in.SetOldRCUFormat(isOldRCUFormat);
 743
 744   // Create a shaper pulse object
 745   AliPHOSPulseGenerator pulse;
 746   TF1 * signalF = new TF1("signal", AliPHOSPulseGenerator::RawResponseFunction, 0, pulse.GetRawFormatTimeMax(), 4);
 747   signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero") ;
 748
 749   Int_t relId[4], absId =0;
 750   Bool_t lowGainFlag = kFALSE ;
 751
 752   TClonesArray * digits = Digits() ;
 753   digits->Clear() ;
 754
 755   Int_t    iBin     = 0;
 756   Int_t    iDigit   = 0 ;
 757   Double_t energyHG = 0. ;
 758   Double_t energyLG = 0. ;
 759   Double_t time     = 0. ;
 760   Int_t    iOldDigit;
 761   Bool_t   seen;
 762   TH1F*    hLowGain  = 0;
 763   TH1F*    hHighGain = 0;
 764
 765   while ( in.Next() ) { // PHOS entries loop
 766
 767     lowGainFlag = in.IsLowGain();
 768
 769     // (Re)create histograms to store samples
 770
 771     if (lowGainFlag) {
 772       if(!hLowGain)
 773         hLowGain = new TH1F("hLowGain","Low gain",in.GetTimeLength(),0,in.GetTimeLength());
 774       else
 775         if(hLowGain->GetNbinsX() != in.GetTimeLength()) {
 776           delete hLowGain;
 777           hLowGain = new TH1F("hLowGain","Low gain",in.GetTimeLength(),0,in.GetTimeLength());
 778         }
 779     } else {
 780       if(!hHighGain)
 781         hHighGain = new TH1F("hHighGain","High gain",in.GetTimeLength(),0,in.GetTimeLength());
 782       else
 783         if(hHighGain->GetNbinsX() != in.GetTimeLength()) {
 784           delete hHighGain;
 785           hHighGain = new TH1F("hHighGain","High gain",in.GetTimeLength(),0,in.GetTimeLength());
 786         }
 787     }
 788
 789     // Fill histograms with samples
 790     if(lowGainFlag)
 791       hLowGain ->SetBinContent(in.GetTimeLength()-iBin-1,in.GetSignal());
 792     else {
 793       hHighGain->SetBinContent(in.GetTimeLength()-iBin-1,in.GetSignal());
 794     }
 795
 796     iBin++;
 797
 798     // Fit the full sample
 799     if(iBin==in.GetTimeLength()) {
 800       iBin=0;
 801
 802       // Temporarily we do not fit the sample graph, but
 803       // take the energy from the graph maximum, and the pedestal from the 0th point
 804       // 30 Aug 2006
 805
 806       //FitRaw(lowGainFlag, hLowGain, hHighGain, signalF, energy, time);
 807
 808       if(lowGainFlag) {
 809         energyLG  = hLowGain ->GetMaximum();     // "digit amplitude"
 810 //      energyLG -= hLowGain ->GetBinContent(0); // "pedestal subtraction"
 811         energyLG *= pulse.GetRawFormatHighLowGainFactor(); // *16
 812         if(AliLog::GetGlobalDebugLevel()>3) {
 813           AliDebug(4,Form("----Printing hLowGain: ----\n")) ;
 814           hLowGain ->Print("all");
 815         }
 816         if(AliLog::GetGlobalDebugLevel()>2)
 817           AliDebug(2,Form("(mod,col,row)=(%d,%d,%d), low gain energy=%f\n\n",
 818                           in.GetModule(),in.GetColumn(),in.GetRow(),energyLG));
 819       }
 820
 821       else {
 822         energyHG  = hHighGain->GetMaximum();     // "digit amplitude"
 823 //      energyHG -= hHighGain->GetBinContent(0); // "pedestal subtraction"
 824         if(AliLog::GetGlobalDebugLevel()>3) {
 825           AliDebug(4,Form("----Printing hHighGain: ----\n"));
 826           hHighGain->Print("all");
 827         }
 828         if(AliLog::GetGlobalDebugLevel()>2)
 829           AliDebug(2,Form("(mod,col,row)=(%d,%d,%d), high gain energy=%f\n\n",
 830                           in.GetModule(),in.GetColumn(),in.GetRow(),energyHG));
 831       }
 832
 833       // Time is not evaluated for the moment (12.01.2007). To be implemented later.
 834       time = -1;
 835
 836       relId[0] = in.GetModule()+1;
 837       relId[1] =                0;
 838       relId[2] = in.GetRow()   +1;
 839       relId[3] = in.GetColumn()+1;
 840       if(!PHOSGeometry()) AliFatal("Couldn't find PHOSGeometry!");
 841       PHOSGeometry()->RelToAbsNumbering(relId, absId);
 842
 843       AliDebug(2,Form("relId=(mod,row,col)=(%d,%d,%d), absId=%d\n\n",
 844                       relId[0],relId[2],relId[3],absId));
 845
 846       // Add low gain digit only if the high gain digit does not exist in the digits array
 847       seen = kFALSE;
 848       if(lowGainFlag) {
 849         for (iOldDigit=iDigit-1; iOldDigit==0; iOldDigit--) {
 850           if ((dynamic_cast<AliPHOSDigit*>(digits->At(iOldDigit)))->GetId() == absId) {
 851             seen = kTRUE;
 852             break;
 853           }
 854         }
 855         if (!seen)
 856           new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyLG,time);
 857         energyLG = 0. ;
 858       }
 859       // Add high gain digit only if it is not saturated;
 860       // replace low gain digit by a high gain one
 861       else {
 862         if (energyHG >= 1023) continue;
 863         for (iOldDigit=iDigit-1; iOldDigit==0; iOldDigit--) {
 864           if ((dynamic_cast<AliPHOSDigit*>(digits->At(iOldDigit)))->GetId() == absId) {
 865             digits->RemoveAt(iOldDigit);
 866             new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyHG,time);
 867             seen = kTRUE;
 868             break;
 869           }
 870         }
 871         if (!seen)
 872           new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyHG,time);
 873         energyHG = 0. ;
 874       }
 875       AliPHOSDigit *dig = dynamic_cast<AliPHOSDigit*>(digits->At(iDigit));
 876       dig->Print();
 877       iDigit++;
 878     }
 879   }
 880
 881   // PHOS entries loop
 882
 883   digits->Sort() ;
 884   digits->Print();
 885
 886   //!!!!for debug!!!
 887   Int_t modMax=-111;
 888   Int_t colMax=-111;
 889   Int_t rowMax=-111;
 890   Float_t eMax=-333;
 891   //!!!for debug!!!
 892
 893   for(iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
 894     AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
 895     if(digit->GetEnergy()>eMax) {
 896       PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
 897       eMax=digit->GetEnergy();
 898       modMax=relId[0];
 899       rowMax=relId[2];
 900       colMax=relId[3];
 901     }
 902   }
 903
 904   AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
 905                   modMax,colMax,rowMax,eMax));
 906
 907   delete signalF ;
 908   delete hHighGain;
 909   delete hLowGain;
 910
 911   return digits->GetEntriesFast() ;
 912 }
 913
 914 //____________________________________________________________________________
 915 Int_t AliPHOSGetter::ReadTreeD() const
 916 {
 917   // Read the Digits
 918
 919   PhosLoader()->CleanDigits() ;
 920   PhosLoader()->LoadDigits("UPDATE") ;
 921   PhosLoader()->LoadDigitizer("UPDATE") ;
 922   return Digits()->GetEntries() ;
 923 }
 924
 925 //____________________________________________________________________________
 926 Int_t AliPHOSGetter::ReadTreeH() const
 927 {
 928   // Read the Hits
 929   PhosLoader()->CleanHits() ;
 930   // gets TreeH from the root file (PHOS.Hit.root)
 931   //if ( !IsLoaded("H") ) {
 932     PhosLoader()->LoadHits("UPDATE") ;
 933   //  SetLoaded("H") ;
 934   //}
 935   return Hits()->GetEntries() ;
 936 }
 937
 938 //____________________________________________________________________________
 939 Int_t AliPHOSGetter::ReadTreeR() const
 940 {
 941   // Read the RecPoints
 942
 943   PhosLoader()->CleanRecPoints() ;
 944   // gets TreeR from the root file (PHOS.RecPoints.root)
 945   //if ( !IsLoaded("R") ) {
 946     PhosLoader()->LoadRecPoints("UPDATE") ;
 947     PhosLoader()->LoadClusterizer("UPDATE") ;
 948     //  SetLoaded("R") ;
 949     //}
 950
 951   return EmcRecPoints()->GetEntries() ;
 952 }
 953
 954 //____________________________________________________________________________
 955 Int_t AliPHOSGetter::ReadTreeT() const
 956 {
 957   // Read the TrackSegments
 958
 959   PhosLoader()->CleanTracks() ;
 960   // gets TreeT from the root file (PHOS.TrackSegments.root)
 961   //if ( !IsLoaded("T") ) {
 962     PhosLoader()->LoadTracks("UPDATE") ;
 963     PhosLoader()->LoadTrackSegmentMaker("UPDATE") ;
 964     //    SetLoaded("T") ;
 965     //}
 966
 967   return TrackSegments()->GetEntries() ;
 968 }
 969 //____________________________________________________________________________
 970 Int_t AliPHOSGetter::ReadTreeP() const
 971 {
 972   // Read the RecParticles
 973
 974   PhosLoader()->CleanRecParticles() ;
 975
 976   // gets TreeT from the root file (PHOS.TrackSegments.root)
 977   //  if ( !IsLoaded("P") ) {
 978     PhosLoader()->LoadRecParticles("UPDATE") ;
 979     PhosLoader()->LoadPID("UPDATE") ;
 980     //  SetLoaded("P") ;
 981     //}
 982
 983   return RecParticles()->GetEntries() ;
 984 }
 985 //____________________________________________________________________________
 986 Int_t AliPHOSGetter::ReadTreeS() const
 987 {
 988   // Read the SDigits
 989
 990   PhosLoader()->CleanSDigits() ;
 991   // gets TreeS from the root file (PHOS.SDigits.root)
 992   //if ( !IsLoaded("S") ) {
 993     PhosLoader()->LoadSDigits("READ") ;
 994     PhosLoader()->LoadSDigitizer("READ") ;
 995     //  SetLoaded("S") ;
 996     //}
 997
 998   return SDigits()->GetEntries() ;
 999 }
1000
1001 //____________________________________________________________________________
1002 Int_t AliPHOSGetter::ReadTreeE(Int_t event)
1003 {
1004   // Read the ESD
1005
1006   // gets esdTree from the root file (AliESDs.root)
1007   if (!fESDFile)
1008     if ( !OpenESDFile() )
1009       return -1 ;
1010
1011   fESDTree = static_cast<TTree*>(fESDFile->Get("esdTree")) ;
1012   fESD = new AliESD;
1013    if (!fESDTree) {
1014
1015      Error("ReadTreeE", "no ESD tree found");
1016      return -1;
1017    }
1018    fESDTree->SetBranchAddress("ESD", &fESD);
1019    fESDTree->GetEvent(event);
1020
1021    return event ;
1022 }
1023
1024 //____________________________________________________________________________
1025 TClonesArray * AliPHOSGetter::SDigits() const
1026 {
1027   // asks the Loader to return the Digits container
1028
1029   TClonesArray * rv = 0 ;
1030
1031   rv = PhosLoader()->SDigits() ;
1032   if (!rv) {
1033     PhosLoader()->MakeSDigitsArray() ;
1034     rv = PhosLoader()->SDigits() ;
1035   }
1036   return rv ;
1037 }
1038
1039 //____________________________________________________________________________
1040 AliPHOSSDigitizer * AliPHOSGetter::SDigitizer()
1041 {
1042   // Returns pointer to the SDigitizer task
1043   AliPHOSSDigitizer * rv ;
1044   rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
1045   if (!rv) {
1046     Event(0, "S") ;
1047     rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
1048   }
1049   return rv ;
1050 }
1051
1052 //____________________________________________________________________________
1053 TParticle * AliPHOSGetter::Secondary(const TParticle* p, Int_t index) const
1054 {
1055   // Return first (index=1) or second (index=2) secondary particle of primary particle p
1056
1057   if(index <= 0)
1058     return 0 ;
1059   if(index > 2)
1060     return 0 ;
1061
1062   if(p) {
1063   Int_t daughterIndex = p->GetDaughter(index-1) ;
1064   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1065   return  rl->GetAliRun()->GetMCApp()->Particle(daughterIndex) ;
1066   }
1067   else
1068     return 0 ;
1069 }
1070
1071 //____________________________________________________________________________
1072 void AliPHOSGetter::Track(Int_t itrack)
1073 {
1074   // Read the first entry of PHOS branch in hit tree gAlice->TreeH()
1075
1076  AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1077
1078   if( !TreeH() ) // load treeH the first time
1079     rl->LoadHits() ;
1080
1081   // first time create the container
1082   TClonesArray * hits = Hits() ;
1083   if ( hits )
1084     hits->Clear() ;
1085
1086   TBranch * phosbranch = dynamic_cast<TBranch*>(TreeH()->GetBranch("PHOS")) ;
1087   phosbranch->SetAddress(&hits) ;
1088   phosbranch->GetEntry(itrack) ;
1089 }
1090
1091 //____________________________________________________________________________
1092 TTree * AliPHOSGetter::TreeD() const
1093 {
1094   // Returns pointer to the Digits Tree
1095   TTree * rv = 0 ;
1096   rv = PhosLoader()->TreeD() ;
1097   if ( !rv ) {
1098     PhosLoader()->MakeTree("D");
1099     rv = PhosLoader()->TreeD() ;
1100   }
1101
1102   return rv ;
1103 }
1104
1105 //____________________________________________________________________________
1106 TTree * AliPHOSGetter::TreeH() const
1107 {
1108   // Returns pointer to the Hits Tree
1109   TTree * rv = 0 ;
1110   rv = PhosLoader()->TreeH() ;
1111   if ( !rv ) {
1112     PhosLoader()->MakeTree("H");
1113     rv = PhosLoader()->TreeH() ;
1114   }
1115
1116   return rv ;
1117 }
1118
1119 //____________________________________________________________________________
1120 TTree * AliPHOSGetter::TreeR() const
1121 {
1122   // Returns pointer to the RecPoints Tree
1123   TTree * rv = 0 ;
1124   rv = PhosLoader()->TreeR() ;
1125   if ( !rv ) {
1126     PhosLoader()->MakeTree("R");
1127     rv = PhosLoader()->TreeR() ;
1128   }
1129
1130   return rv ;
1131 }
1132
1133 //____________________________________________________________________________
1134 TTree * AliPHOSGetter::TreeT() const
1135 {
1136   // Returns pointer to the TrackSegments Tree
1137   TTree * rv = 0 ;
1138   rv = PhosLoader()->TreeT() ;
1139   if ( !rv ) {
1140     PhosLoader()->MakeTree("T");
1141     rv = PhosLoader()->TreeT() ;
1142   }
1143
1144   return rv ;
1145 }
1146 //____________________________________________________________________________
1147 TTree * AliPHOSGetter::TreeP() const
1148 {
1149   // Returns pointer to the RecParticles  Tree
1150   TTree * rv = 0 ;
1151   rv = PhosLoader()->TreeP() ;
1152   if ( !rv ) {
1153     PhosLoader()->MakeTree("P");
1154     rv = PhosLoader()->TreeP() ;
1155   }
1156
1157   return rv ;
1158 }
1159
1160 //____________________________________________________________________________
1161 TTree * AliPHOSGetter::TreeS() const
1162 {
1163  // Returns pointer to the SDigits Tree
1164   TTree * rv = 0 ;
1165   rv = PhosLoader()->TreeS() ;
1166   if ( !rv ) {
1167     PhosLoader()->MakeTree("S");
1168     rv = PhosLoader()->TreeS() ;
1169   }
1170
1171   return rv ;
1172 }
1173
1174 //____________________________________________________________________________
1175 Bool_t AliPHOSGetter::VersionExists(TString & opt) const
1176 {
1177   // checks if the version with the present name already exists in the same directory
1178
1179   Bool_t rv = kFALSE ;
1180
1181   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1182   TString version( rl->GetEventFolder()->GetName() ) ;
1183
1184   opt.ToLower() ;
1185
1186   if ( opt == "sdigits") {
1187     // add the version name to the root file name
1188     TString fileName( PhosLoader()->GetSDigitsFileName() ) ;
1189     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1190       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1191     if ( !(gSystem->AccessPathName(fileName)) ) {
1192       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1193       rv = kTRUE ;
1194     }
1195     PhosLoader()->SetSDigitsFileName(fileName) ;
1196   }
1197
1198   if ( opt == "digits") {
1199     // add the version name to the root file name
1200     TString fileName( PhosLoader()->GetDigitsFileName() ) ;
1201     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1202       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1203     if ( !(gSystem->AccessPathName(fileName)) ) {
1204       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1205       rv = kTRUE ;
1206     }
1207   }
1208
1209   return rv ;
1210
1211 }
1212
1213 //____________________________________________________________________________
1214 UShort_t AliPHOSGetter::EventPattern(void) const
1215 {
1216   // Return the pattern (trigger bit register) of the beam-test event
1217   if(fBTE)
1218     return fBTE->GetPattern() ;
1219   else
1220     return 0 ;
1221 }
1222 //____________________________________________________________________________
1223 Float_t AliPHOSGetter::BeamEnergy(void) const
1224 {
1225   // Return the beam energy of the beam-test event
1226   if(fBTE)
1227     return fBTE->GetBeamEnergy() ;
1228   else
1229     return 0 ;
1230 }
1231 //____________________________________________________________________________
1232
1233 AliPHOSCalibData* AliPHOSGetter::CalibData()
1234 {
1235   // Check if the instance of AliPHOSCalibData exists, and return it
1236
1237   return fgCalibData;
1238 }