Adding statistical function (Marian)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSSDigitizer.cxx
index b0ca76b..bcffe97 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+
 /* $Id$ */
 
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.49  2006/05/10 06:42:53  kharlov
+ * Remove redundant loop over primaries
+ *
+ * Revision 1.48  2006/04/22 10:30:17  hristov
+ * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
+ *
+ * Revision 1.47  2005/05/28 14:19:04  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //_________________________________________________________________________
 // This is a TTask that makes SDigits out of Hits
+// The name of the TTask is also the title of the branch that will contain 
+// the created SDigits
+// The title of the TTAsk is the name of the file that contains the hits from
+// which the SDigits are created
 // A Summable Digits is the sum of all hits originating 
 // from one primary in one active cell
 // A threshold for assignment of the primary to SDigit is applied 
 
 
 // --- ROOT system ---
-#include "TFile.h"
-#include "TTask.h"
-#include "TTree.h"
-#include "TSystem.h"
-#include "TROOT.h"
-#include "TFolder.h"
 #include "TBenchmark.h"
+               //#include "TObjectTable.h"
+
 // --- Standard library ---
-#include <iomanip.h>
 
 // --- AliRoot header files ---
-#include "AliRun.h"
+#include "AliLog.h"
+#include "AliPHOSGeometry.h" 
 #include "AliPHOSDigit.h"
+#include "AliPHOSGetter.h"
 #include "AliPHOSHit.h"
 #include "AliPHOSSDigitizer.h"
-
+//#include "AliMemoryWatcher.h"
 
 ClassImp(AliPHOSSDigitizer)
 
            
 //____________________________________________________________________________ 
-  AliPHOSSDigitizer::AliPHOSSDigitizer():TTask("AliPHOSSDigitizer","") 
+AliPHOSSDigitizer::AliPHOSSDigitizer() : 
+  TTask("",""),
+  fA(0.f), fB(0.f),
+  fPrimThreshold(0.f),
+  fDefaultInit(kTRUE),
+  fEventFolderName(""),
+  fInit(kFALSE),
+  fSDigitsInRun(0),
+  fFirstEvent(0),
+  fLastEvent(0)
 {
   // ctor
-  fA             = 0;
-  fB             = 10000000.;
-  fPrimThreshold = 0.01 ;
-  fNevents       = 0 ;      
-  fSDigits       = 0 ;
-  fHits          = 0 ;
-  fIsInitialized = kFALSE ;
-
 }
 
 //____________________________________________________________________________ 
-AliPHOSSDigitizer::AliPHOSSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliPHOSSDigitizer","")
+AliPHOSSDigitizer::AliPHOSSDigitizer(const char * alirunFileName, 
+                                    const char * eventFolderName):
+  TTask("PHOS"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
+  fA(0.f), fB(0.f),
+  fPrimThreshold(0.f),
+  fDefaultInit(kFALSE),
+  fEventFolderName(eventFolderName),
+  fInit(kFALSE),
+  fSDigitsInRun(0),
+  fFirstEvent(0),
+  fLastEvent(0)
 {
   // ctor
-  fA             = 0;
-  fB             = 10000000.;
-  fPrimThreshold = 0.01 ;
-  fNevents       = 0 ;      
-  fSDigitsTitle  = sDigitsTitle ;
-  fHeadersFile   = headerFile ;
-  fIsInitialized = kFALSE ;
-
+  InitParameters() ; 
   Init();
+  fDefaultInit = kFALSE ; 
 }
 
 //____________________________________________________________________________ 
-AliPHOSSDigitizer::~AliPHOSSDigitizer()
-{
-  // dtor
-  if(fSDigits)
-    delete fSDigits ;
-  if(fHits)
-    delete fHits ;
+AliPHOSSDigitizer::~AliPHOSSDigitizer() {
+  //dtor
+  //  AliPHOSGetter * gime =
+  //  AliPHOSGetter::Instance(GetTitle(),fEventFolderName.Data());  
+  AliPHOSGetter * gime =
+    AliPHOSGetter::Instance();  
+  gime->PhosLoader()->CleanSDigitizer();
 }
 //____________________________________________________________________________ 
 void AliPHOSSDigitizer::Init()
 {
-  // Initialization: open root-file, allocate arrays for hits and sdigits,
-  // attach task SDigitizer to the list of PHOS tasks
-  // 
-  // Initialization can not be done in the default constructor
-
-  if(!fIsInitialized){
+  // Uses the getter to access the required files
+  
+  fInit = kTRUE ; 
+  
+  AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());  
+  if ( gime == 0 ) {
+    Fatal("Init" ,"Could not obtain the Getter object for file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;  
+    return ;
+  } 
+  
+  TString opt("SDigits") ; 
+  if(gime->VersionExists(opt) ) { 
+    Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
+    fInit = kFALSE ; 
+  }
 
-    if(fHeadersFile.IsNull())
-      fHeadersFile="galice.root" ;
+  gime->PostSDigitizer(this);
+  gime->PhosLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
+}
 
-    TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
-    
-    //if file was not opened yet, read gAlice
-    if(file == 0){
-      if(fHeadersFile.Contains("rfio"))
-       file =  TFile::Open(fHeadersFile,"update") ;
-      else
-       file = new TFile(fHeadersFile.Data(),"update") ;
-      gAlice = (AliRun *) file->Get("gAlice") ;
-    }
-    
-    fHits    = new TClonesArray("AliPHOSHit",1000);
-    fSDigits = new TClonesArray("AliPHOSDigit",1000);
-    
-    //add Task to //YSAlice/tasks/(S)Diditizer/PHOS
-    TFolder * alice  = (TFolder*)gROOT->GetListOfBrowsables()->FindObject("YSAlice") ; 
-    TTask * aliceSD  = (TTask*)alice->FindObject("tasks/(S)Digitizer") ; 
-    TTask * phosSD   = (TTask*)aliceSD->GetListOfTasks()->FindObject("PHOS") ;
-    phosSD->Add(this) ; 
-    
-    fIsInitialized = kTRUE ;
-  }
+//____________________________________________________________________________ 
+void AliPHOSSDigitizer::InitParameters()
+{ 
+  // initializes the parameters for digitization
+  fA             = 0;
+  fB             = 10000000.;
+  fPrimThreshold = 0.01 ;
+  fSDigitsInRun  = 0 ;
 }
+
 //____________________________________________________________________________
 void AliPHOSSDigitizer::Exec(Option_t *option) 
 { 
-  // Collects all hits in the same active volume into digit
+  // Steering method to produce summable digits for events
+  // in the range from fFirstEvent to fLastEvent.
+  // This range is optionally set by SetEventRange().
+  // if fLastEvent=-1 (by default), then process events until the end.
+  //
+  // Summable digit is a sum of all hits in the same active
+  // volume into digit
   
-  if(!fIsInitialized)
-    Init() ;
+  if (strstr(option, "print") ) {
+    Print() ; 
+    return ; 
+  }
 
   if(strstr(option,"tim"))
     gBenchmark->Start("PHOSSDigitizer");
   
-  fNevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
+
+  //switch off reloading of this task while getting event
+  if (!fInit) { // to prevent overwrite existing file
+    AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
+    return ;
+  }
+
+  if (fLastEvent == -1) 
+    fLastEvent = gime->MaxEvent() - 1 ;
+  else 
+    fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // only ine event at the time 
+  Int_t nEvents   = fLastEvent - fFirstEvent + 1;
   
-  Int_t ievent ;
-  for(ievent = 0; ievent < fNevents; ievent++){
-    gAlice->GetEvent(ievent) ;
-    gAlice->SetEvent(ievent) ;
-    
-    if(gAlice->TreeH()==0){
-      cout << "AliPHOSSDigitizer: There is no Hit Tree" << endl;
-      return ;
-    }
-    
-    //set address of the hits 
-    TBranch * branch = gAlice->TreeH()->GetBranch("PHOS");
-    if (branch) 
-      branch->SetAddress(&fHits);
-    else{
-      cout << "ERROR in AliPHOSSDigitizer: "<< endl ;
-      cout << "      no branch PHOS in TreeH"<< endl ;
-      cout << "      do nothing " << endl ;
-      return ;
-    }
-    
-    fSDigits->Clear();
+  Int_t ievent, i;
+
+  //AliMemoryWatcher memwatcher;
+
+  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {    
+    gime->Event(ievent,"H") ;   
+    TTree * treeS = gime->TreeS(); 
+    TClonesArray * hits = gime->Hits() ;
+    TClonesArray * sdigits = gime->SDigits() ;
+    sdigits->Clear();
     Int_t nSdigits = 0 ;
-    
-    
-    //Now made SDigits from hits, for PHOS it is the same, so just copy    
-    Int_t itrack ;
-    for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){
-      
-      //=========== Get the PHOS branch from Hits Tree for the Primary track itrack
-      branch->GetEntry(itrack,0);
+    //Now make SDigits from hits, for PHOS it is the same, so just copy    
+    for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
+      AliPHOSHit * hit = dynamic_cast<AliPHOSHit *>(hits->At(i)) ;
+      // Assign primary number only if contribution is significant
       
-      Int_t i;
-      for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
-       AliPHOSHit * hit = (AliPHOSHit*)fHits->At(i) ;
-
-       // Assign primary number only if contribution is significant
-       if( hit->GetEnergy() > fPrimThreshold)
-         new((*fSDigits)[nSdigits]) AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-       else
-         new((*fSDigits)[nSdigits]) AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-       
-       nSdigits++ ;  
-       
-      } 
+      if( hit->GetEnergy() > fPrimThreshold)
+       new((*sdigits)[nSdigits]) AliPHOSDigit(hit->GetPrimary(),hit->GetId(),
+                                              hit->GetEnergy() ,hit->GetTime()) ;
+      else
+       new((*sdigits)[nSdigits]) AliPHOSDigit(-1               ,hit->GetId(), 
+                                              hit->GetEnergy() ,hit->GetTime()) ;
+      nSdigits++ ;     
       
-    } // loop over tracks
-    
-    fSDigits->Sort() ;
-    
-    nSdigits = fSDigits->GetEntriesFast() ;
-    fSDigits->Expand(nSdigits) ;
-    
-    Int_t i ;
+    }
+    sdigits->Sort() ;
+
+    nSdigits = sdigits->GetEntriesFast() ;
+
+    fSDigitsInRun += nSdigits ;  
+    sdigits->Expand(nSdigits) ;
+
     for (i = 0 ; i < nSdigits ; i++) { 
-      AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ; 
+      AliPHOSDigit * digit = dynamic_cast<AliPHOSDigit *>(sdigits->At(i)) ; 
       digit->SetIndexInList(i) ;     
     }
 
-    if(gAlice->TreeS() == 0)
-      gAlice->MakeTree("S") ;
-    
-    //check, if this branch already exits?
-    TBranch * sdigitsBranch = 0;
-    TBranch * sdigitizerBranch = 0;
-    
-    TObjArray * branches = gAlice->TreeS()->GetListOfBranches() ;
-    Int_t ibranch;
-    Bool_t phosNotFound = kTRUE ;
-    Bool_t sdigitizerNotFound = kTRUE ;
-    
-    for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
-      
-      if(phosNotFound){
-       sdigitsBranch=(TBranch *) branches->At(ibranch) ;
-       if( (strcmp("PHOS",sdigitsBranch->GetName())==0 ) &&
-           (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) )
-         phosNotFound = kFALSE ;
-      }
-      if(sdigitizerNotFound){
-       sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
-       if( (strcmp(sdigitizerBranch->GetName(),"AliPHOSSDigitizer") == 0)&&
-           (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) )
-         sdigitizerNotFound = kFALSE ;
-      }
-    }
+    //Now write SDigits
 
-    if(!(sdigitizerNotFound && phosNotFound)){
-      cout << "AliPHOSSdigitizer error:" << endl ;
-      cout << "Can not overwrite existing branches: do not write" << endl ;
-      return ;
-    }
-    
-    //Make (if necessary) branches    
-    char * file =0;
-    if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
-      file = new char[strlen(gAlice->GetBaseFile())+20] ;
-      sprintf(file,"%s/PHOS.SDigits.root",gAlice->GetBaseFile()) ;
-    }
-    
-    TDirectory *cwd = gDirectory;
     
     //First list of sdigits
-    Int_t bufferSize = 32000 ;    
-    sdigitsBranch = gAlice->TreeS()->Branch("PHOS",&fSDigits,bufferSize);
-    sdigitsBranch->SetTitle(fSDigitsTitle.Data());
-    if (file) {
-      sdigitsBranch->SetFile(file);
-      TIter next( sdigitsBranch->GetListOfBranches());
-      TBranch * subbr;
-      while ((subbr=(TBranch*)next())) {
-       subbr->SetFile(file);
-      }   
-      cwd->cd();
-    } 
-      
-    //second - SDigitizer
-    Int_t splitlevel = 0 ;
-    AliPHOSSDigitizer * sd = this ;
-    sdigitizerBranch = gAlice->TreeS()->Branch("AliPHOSSDigitizer","AliPHOSSDigitizer",
-                                              &sd,bufferSize,splitlevel); 
-    sdigitizerBranch->SetTitle(fSDigitsTitle.Data());
-    if (file) {
-      sdigitizerBranch->SetFile(file);
-      TIter next( sdigitizerBranch->GetListOfBranches());
-      TBranch * subbr ;
-      while ((subbr=(TBranch*)next())) {
-       subbr->SetFile(file);
-      }   
-      cwd->cd();
-      delete file;
-    }
+
+    Int_t bufferSize = 32000 ;
+    TBranch * sdigitsBranch = treeS->Branch("PHOS",&sdigits,bufferSize);
 
     sdigitsBranch->Fill() ;
-    sdigitizerBranch->Fill() ;
-    gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
-    
+
+    gime->WriteSDigits("OVERWRITE");
+
+    //Next - SDigitizer
+
+    gime->WriteSDigitizer("OVERWRITE");
+    //gObjectTable->Print() ; 
+  
     if(strstr(option,"deb"))
       PrintSDigits(option) ;
-    
-  }
+      
+    //memwatcher.Watch(ievent); 
+  }// event loop
+  
+  Unload();
+
+  //  gime->PhosLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
   
   if(strstr(option,"tim")){
     gBenchmark->Stop("PHOSSDigitizer");
-    cout << "AliPHOSSDigitizer:" << endl ;
-    cout << "   took " << gBenchmark->GetCpuTime("PHOSSDigitizer") << " seconds for SDigitizing " 
-        <<  gBenchmark->GetCpuTime("PHOSSDigitizer")/fNevents << " seconds per event " << endl ;
-    cout << endl ;
+    Info("Exec","   took %f seconds for SDigitizing  %f seconds per event",
+        gBenchmark->GetCpuTime("PHOSSDigitizer"), gBenchmark->GetCpuTime("PHOSSDigitizer")/nEvents) ;
   }
   
-  
-}
-//__________________________________________________________________
-void AliPHOSSDigitizer::SetSDigitsBranch(const char * title )
-{
-  // Setting title to branch SDigits 
-  if(!fSDigitsTitle.IsNull())
-    cout << "AliPHOSSdigitizer: changing SDigits file from " <<fSDigitsTitle.Data() << " to " << title << endl ;
-  fSDigitsTitle=title ;
+  //TFile f("out.root","RECREATE");
+  //memwatcher.WriteToFile(); 
+  //f.Close();
 }
+
 //__________________________________________________________________
-void AliPHOSSDigitizer::Print(Option_t* option)const
+void AliPHOSSDigitizer::Print(const Option_t *)const
 {
   // Prints parameters of SDigitizer
-  cout << "------------------- "<< GetName() << " -------------" << endl ;
-  cout << "   Writing SDigitis to branch with title  " << fSDigitsTitle.Data() << endl ;
-  cout << "   with digitization parameters  A = " << fA << endl ;
-  cout << "                                 B = " << fB << endl ;
-  cout << "   Threshold for Primary assignment= " << fPrimThreshold << endl ; 
-  cout << "---------------------------------------------------"<<endl ;
+  Info("Print", "\n------------------- %s -------------", GetName() ) ; 
+  printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
+  printf("   with digitization parameters  A = %f\n", fA) ; 
+  printf("                                 B = %f\n", fB) ;
+  printf("   Threshold for Primary assignment= %f\n", fPrimThreshold)  ; 
+  printf("---------------------------------------------------\n") ;
   
 }
+
 //__________________________________________________________________
 Bool_t AliPHOSSDigitizer::operator==( AliPHOSSDigitizer const &sd )const
 {
@@ -336,32 +292,67 @@ Bool_t AliPHOSSDigitizer::operator==( AliPHOSSDigitizer const &sd )const
   else
     return kFALSE ;
 }
+
 //__________________________________________________________________
 void AliPHOSSDigitizer::PrintSDigits(Option_t * option)
 {
   // Prints list of digits produced in the current pass of AliPHOSDigitizer
+
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  const TClonesArray * sdigits = gime->SDigits() ;
   
-  cout << "AliPHOSSDigitizer: " << endl ;
-  cout << "       Number of entries in SDigits list  " << fSDigits->GetEntriesFast() << endl ;
-  cout << endl ;
-  
-  if(strstr(option,"all")){// print all digits
+  Info( "\nPrintSDigits", "event # %d %d sdigits", gAlice->GetEvNumber(), sdigits->GetEntriesFast() ) ; 
+
+  if(strstr(option,"all")||strstr(option,"EMC")){
     
     //loop over digits
     AliPHOSDigit * digit;
-    cout << "SDigit Id " << " Amplitude " <<  " Index "  <<  " Nprim " << " Primaries list " <<  endl;    
+    printf("\nEMC sdigits\n") ; 
+    Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
     Int_t index ;
-    for (index = 0 ; index < fSDigits->GetEntries() ; index++) {
-      digit = (AliPHOSDigit * )  fSDigits->At(index) ;
-      cout << setw(8)  <<  digit->GetId() << " "  <<   setw(3)  <<  digit->GetAmp() <<   "  "  
-          << setw(6)  <<  digit->GetIndexInList() << "  "   
-          << setw(5)  <<  digit->GetNprimary() <<"  ";
-      
+    for (index = 0 ; (index < sdigits->GetEntriesFast()) && 
+        ((dynamic_cast<AliPHOSDigit *> (sdigits->At(index)))->GetId() <= maxEmc) ; index++) {
+      digit = dynamic_cast<AliPHOSDigit *>( sdigits->At(index) ) ;
+      //  if(digit->GetNprimary() == 0) 
+      //       continue;
+//       printf("%6d  %8d    %6.5e %4d      %2d :\n", // YVK
+      printf("%6d  %.4f    %6.5e %4d      %2d :\n",
+             digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
       Int_t iprimary;
-      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
-       cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "  ";
-      cout << endl;     
-    }
+      for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
+       printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
+      }  
+    }    
+  }
+
+  if(strstr(option,"all")||strstr(option,"CPV")){
     
+    //loop over CPV digits
+    AliPHOSDigit * digit;
+    printf("\nCPV sdigits\n") ; 
+    Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
+    Int_t index ;
+    for (index = 0 ; index < sdigits->GetEntriesFast(); index++) {
+      digit = dynamic_cast<AliPHOSDigit *>( sdigits->At(index) ) ;
+      if(digit->GetId() > maxEmc){
+       printf("\n%6d  %8d    %4d      %2d :",
+               digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;  
+       Int_t iprimary;
+       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
+         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
+       }
+      }    
+    }
   }
 }
+
+//____________________________________________________________________________ 
+void AliPHOSSDigitizer::Unload() const
+{
+  // Unloads the objects from the folder
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  AliPHOSLoader * loader = gime->PhosLoader() ; 
+  loader->UnloadHits() ; 
+  loader->UnloadSDigits() ; 
+}