Adding statistical function (Marian)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSSDigitizer.cxx
index 98dff99..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 
 
 
 // --- 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() : 
+  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 ;
-  fSDigitsInRun  = 0 ; 
 }
 
 //____________________________________________________________________________ 
-AliPHOSSDigitizer::AliPHOSSDigitizer(const char * headerFile, const char * sDigitsTitle):TTask(sDigitsTitle, headerFile)
+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 ;
-  fSDigitsInRun  = 0 ; 
+  InitParameters() ; 
   Init();
+  fDefaultInit = kFALSE ; 
 }
 
-
+//____________________________________________________________________________ 
+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
-  //============================================================= YS
-  //  The initialisation is now done by AliPHOSGetter
+  // Uses the getter to access the required files
   
-  if( strcmp(GetTitle(), "") == 0 )
-    SetTitle("galice.root") ;
+  fInit = kTRUE ; 
   
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), GetName()) ;     
+  AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());  
   if ( gime == 0 ) {
-    cerr << "ERROR: AliPHOSSDigitizer::Init -> Could not obtain the Getter object !" << endl ; 
+    Fatal("Init" ,"Could not obtain the Getter object for file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;  
     return ;
   } 
   
-  //  //add Task to //YSAlice/tasks/SDiditizer/PHOS
-  //  TTask * aliceSD  = (TTask*)gROOT->FindObjectAny("YSAlice/tasks/SDigitizer") ; 
-  //  TTask * phosSD   = (TTask*)aliceSD->GetListOfTasks()->FindObject("PHOS") ;
-  //  if ( ! phosSD->GetListOfTasks()->FindObject(GetName()) ) 
-  //   phosSD->Add(this) ;
-
-  gime->PostSDigits( GetName(), GetTitle() ) ; 
-
-  TString sdname(GetName() );
-  sdname.Append(":") ;
-  sdname.Append(GetTitle() ) ;
-  SetName(sdname) ;
-  gime->PostSDigitizer(this) ;
+  TString opt("SDigits") ; 
+  if(gime->VersionExists(opt) ) { 
+    Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
+    fInit = kFALSE ; 
+  }
+
+  gime->PostSDigitizer(this);
+  gime->PhosLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
  
-  // create a folder on the white board //YSAlice/WhiteBoard/SDigits/PHOS/headerFile/sdigitsTitle
-    
+}
+
+//____________________________________________________________________________ 
+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
-  if( strcmp(GetName(), "") == 0 )
-    Init() ;
+  // 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 (strstr(option, "print") ) {
-    Print("") ; 
+    Print() ; 
     return ; 
   }
 
   if(strstr(option,"tim"))
     gBenchmark->Start("PHOSSDigitizer");
+  
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 
-  //Check, if this branch already exits
-  gAlice->GetEvent(0) ;
-  if(gAlice->TreeS() ) {
-    TObjArray * lob = (TObjArray*)gAlice->TreeS()->GetListOfBranches() ;
-    TIter next(lob) ; 
-    TBranch * branch = 0 ;  
-    Bool_t phosfound = kFALSE, sdigitizerfound = kFALSE ; 
-    
-    while ( (branch = (TBranch*)next()) && (!phosfound || !sdigitizerfound) ) {
-      if ( (strcmp(branch->GetName(), "PHOS")==0) && (strcmp(branch->GetTitle(), GetName())==0) ) 
-       phosfound = kTRUE ;
-      
-      else if ( (strcmp(branch->GetName(), "AliPHOSSDigitizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) ) 
-       sdigitizerfound = kTRUE ; 
-    }
-    
-    if ( phosfound || sdigitizerfound ) {
-      cerr << "WARNING: AliPHOSSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName() 
-          << " already exits" << endl ;
-      return ; 
-    }   
-  }  
-
-  TString sdname(GetName()) ;
-  sdname.Remove(sdname.Index(GetTitle())-1) ;
-    
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
+  //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 ;
+  }
 
-  Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
+  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 < nevents; ievent++){
-    gime->Event(ievent,"H") ;
+  Int_t ievent, i;
+
+  //AliMemoryWatcher memwatcher;
 
-    TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
+  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 make 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
-      gime->Track(itrack) ;
-      TClonesArray * hits = gime->Hits() ;
-      Int_t i;
-      for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
-       AliPHOSHit * hit = (AliPHOSHit *) hits->At(i) ;
-       // Assign primary number only if contribution is significant
-
-       if( hit->GetEnergy() > fPrimThreshold)
-         new((*sdigits)[nSdigits]) AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-       else
-         new((*sdigits)[nSdigits]) AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-       
-       nSdigits++ ;  
-       
-      } 
-
-    } // loop over tracks
-    
+    for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
+      AliPHOSHit * hit = dynamic_cast<AliPHOSHit *>(hits->At(i)) ;
+      // Assign primary number only if contribution is significant
+      
+      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++ ;     
+      
+    }
     sdigits->Sort() ;
-    
+
     nSdigits = sdigits->GetEntriesFast() ;
 
+    fSDigitsInRun += nSdigits ;  
     sdigits->Expand(nSdigits) ;
-    Int_t i ;
+
     for (i = 0 ; i < nSdigits ; i++) { 
-      AliPHOSDigit * digit = (AliPHOSDigit *) sdigits->At(i) ; 
+      AliPHOSDigit * digit = dynamic_cast<AliPHOSDigit *>(sdigits->At(i)) ; 
       digit->SetIndexInList(i) ;     
     }
 
-    if(gAlice->TreeS() == 0)
-      gAlice->MakeTree("S") ;
-    
-    //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;
+    //Now write SDigits
+
     
     //First list of sdigits
-    Int_t bufferSize = 32000 ;    
-    TBranch * sdigitsBranch = gAlice->TreeS()->Branch("PHOS",&sdigits,bufferSize);
-    sdigitsBranch->SetTitle(sdname);
-    if (file) {
-      sdigitsBranch->SetFile(file);
-      TIter next( sdigitsBranch->GetListOfBranches());
-      TBranch * subbr;
-      while ((subbr=(TBranch*)next())) {
-       subbr->SetFile(file);
-      }   
-      cwd->cd();
-    } 
-      
-    //Next - SDigitizer
-    Int_t splitlevel = 0 ;
-    AliPHOSSDigitizer * sd = this ;
-    TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliPHOSSDigitizer","AliPHOSSDigitizer",
-                                              &sd,bufferSize,splitlevel); 
-    sdigitizerBranch->SetTitle(sdname);
-    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")/nevents << " seconds per event " << endl ;
-    cout << endl ;
-  }
-  
-  
-}
-//__________________________________________________________________
-void AliPHOSSDigitizer::SetSDigitsBranch(const char * title )
-{
-  // Setting title to branch SDigits 
-
-  TString stitle(title) ;
-
-  // check if branch with title already exists
-  TBranch * sdigitsBranch    =  (TBranch*)gAlice->TreeS()->GetListOfBranches()->FindObject("PHOS") ; 
-  TBranch * sdigitizerBranch =  (TBranch*)gAlice->TreeS()->GetListOfBranches()->FindObject("AliPHOSSDigitizer") ;
-  const char * sdigitsTitle    = sdigitsBranch ->GetTitle() ;  
-  const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
-  if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
-    cerr << "ERROR: AliPHOSSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
-    return ;
+    Info("Exec","   took %f seconds for SDigitizing  %f seconds per event",
+        gBenchmark->GetCpuTime("PHOSSDigitizer"), gBenchmark->GetCpuTime("PHOSSDigitizer")/nEvents) ;
   }
   
-  cout << "AliPHOSSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
-
-  SetName(title) ; 
-    
-  // Post to the WhiteBoard
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  gime->PostSDigits( title, GetTitle()) ; 
+  //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  " << GetName() << 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
 {
@@ -329,46 +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::GetInstance() ; 
-  TString sdname(GetName()) ;
-  sdname.Remove(sdname.Index(GetTitle())-1) ;
-  TClonesArray * sdigits = gime->SDigits(sdname.Data()) ; 
-
-  cout << "AliPHOSSDigitizer: event "<<gAlice->GetEvNumber() << endl ;
-  cout << "       Number of entries in SDigits list  " << sdigits->GetEntriesFast() << endl ;
-  cout << endl ;
 
-  fSDigitsInRun +=  sdigits->GetEntriesFast() ; 
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  const TClonesArray * sdigits = gime->SDigits() ;
   
-  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 < sdigits->GetEntries() ; index++) {
-      digit = (AliPHOSDigit * )  sdigits->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::UseHitsFrom(const char * filename)
+void AliPHOSSDigitizer::Unload() const
 {
-  SetTitle(filename) ; 
-  Init() ; 
+  // Unloads the objects from the folder
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  AliPHOSLoader * loader = gime->PhosLoader() ; 
+  loader->UnloadHits() ; 
+  loader->UnloadSDigits() ; 
 }