]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSSDigitizer.cxx
Adding statistical function (Marian)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSSDigitizer.cxx
index 5cd484e927f58955cd3600625acdde887de92706..bcffe972ec534f07e02cc8ef7e57cb171eee9711 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 "TGeometry.h"
+               //#include "TObjectTable.h"
 
 // --- Standard library ---
-#include <iomanip.h>
 
 // --- AliRoot header files ---
-#include "AliRun.h"
-#include "AliHeader.h"
+#include "AliLog.h"
+#include "AliPHOSGeometry.h" 
 #include "AliPHOSDigit.h"
-#include "AliPHOSGeometry.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 ;
   } 
   
-  gime->PostSDigits( GetName(), GetTitle() ) ; 
+  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);
+}
 
-  TString sdname(GetName() );
-  sdname.Append(":") ;
-  sdname.Append(GetTitle() ) ;
-  SetName(sdname) ;
-  gime->PostSDigitizer(this) ;
+//____________________________________________________________________________ 
+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) ;
+  //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 ;
+  }
 
-  TString sdname(GetName()) ;
-  sdname.Remove(sdname.Index(GetTitle())-1) ;
+  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;
   
-  if(gAlice->TreeS() ) {
-    TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
-    TIter next(lob) ; 
-    TBranch * branch = 0 ;  
-    Bool_t phosfound = kFALSE, sdigitizerfound = kFALSE ; 
-    
-    while ( (branch = (static_cast<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 ; 
-    }   
-  }  
+  Int_t ievent, i;
 
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
+  //AliMemoryWatcher memwatcher;
 
-  Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
-  
-  Int_t ievent ;
-  for(ievent = 0; ievent < nevents; ievent++){
-    gime->Event(ievent,"H") ;
-    const TClonesArray * hits = gime->Hits() ;
-    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 Nprim =  (Int_t)  (gAlice->TreeH())->GetEntries(); 
-    // Attention Nprim is the number of primaries tracked by Geant and this number could be different to the number of Primaries in TreeK;
-    Int_t iprim;
-    for (iprim=0; iprim<Nprim; iprim++) { 
-      //=========== Get the PHOS branch from Hits Tree for the Primary iprim
-      gime->Track(iprim) ;
-      Int_t i;
-      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(),
-                                                 Digitize(hit->GetEnergy()), hit->GetTime()) ;
-       else
-         new((*sdigits)[nSdigits]) AliPHOSDigit( -1              , hit->GetId(), 
-                                                  Digitize(hit->GetEnergy()), hit->GetTime()) ;
-       nSdigits++ ;    
-       
-      }
-    } // loop over iprim
-    
+    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 = dynamic_cast<AliPHOSDigit *>(sdigits->At(i)) ; 
       digit->SetIndexInList(i) ;     
     }
 
-    if(gAlice->TreeS() == 0)
-       gAlice->MakeTree("S") ;
-        
+    //Now write SDigits
+
+    
     //First list of sdigits
-    Int_t bufferSize = 32000 ;    
-    TBranch * sdigitsBranch = gAlice->TreeS()->Branch("PHOS",&sdigits,bufferSize);
-    sdigitsBranch->SetTitle(sdname);
-      
+
+    Int_t bufferSize = 32000 ;
+    TBranch * sdigitsBranch = treeS->Branch("PHOS",&sdigits,bufferSize);
+
+    sdigitsBranch->Fill() ;
+
+    gime->WriteSDigits("OVERWRITE");
+
     //Next - SDigitizer
-    Int_t splitlevel = 0 ;
-    AliPHOSSDigitizer * sd = this ;
-    TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliPHOSSDigitizer","AliPHOSSDigitizer",
-                                              &sd,bufferSize,splitlevel); 
-    sdigitizerBranch->SetTitle(sdname);
-    sdigitsBranch->Fill() ; 
-
-    gAlice->TreeS()->AutoSave() ;
-    
+
+    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    = 
-    static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("PHOS")) ; 
-  TBranch * sdigitizerBranch =  
-    static_cast<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()) ; 
-}
-
-//__________________________________________________________________
-void AliPHOSSDigitizer::SetSplitFile(const TString splitFileName) const
-{
-  // Diverts the SDigits in a file separate from the hits file
-  
-  TDirectory * cwd = gDirectory ;
-  TFile * splitFile = gAlice->InitTreeFile("S",splitFileName.Data());
-  splitFile->cd() ; 
-  gAlice->Write();
-  
-  TTree *treeE  = gAlice->TreeE();
-  if (!treeE) {
-    cerr<<"No TreeE found "<<endl;
-    abort() ;
-  }      
-  
-  // copy TreeE
-  AliHeader *header = new AliHeader();
-  treeE->SetBranchAddress("Header", &header);
-  treeE->SetBranchStatus("*",1);
-  TTree *treeENew =  treeE->CloneTree();
-  treeENew->Write();
-  
-  // copy AliceGeom
-  TGeometry *AliceGeom = static_cast<TGeometry*>(cwd->Get("AliceGeom"));
-  if (!AliceGeom) {
-    cerr<<"AliceGeom was not found in the input file "<<endl;
-    abort() ;
-  }
-  AliceGeom->Write();
-  cwd->cd() ; 
-  gAlice->MakeTree("S",splitFile);
-  cout << "INFO: AliPHOSSDigitizer::SetSPlitMode -> SDigits will be stored in " << splitFileName.Data() << endl ; 
+  //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 SDigits 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
 {
@@ -334,73 +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) ;
-  const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ; 
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  const TClonesArray * sdigits = gime->SDigits() ;
+  
+  Info( "\nPrintSDigits", "event # %d %d sdigits", gAlice->GetEvNumber(), sdigits->GetEntriesFast() ) ; 
 
-  cout << "AliPHOSSDigitiser: event " << gAlice->GetEvNumber() << endl ;
-  cout << "      Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
-  cout << endl ;
   if(strstr(option,"all")||strstr(option,"EMC")){
     
     //loop over digits
     AliPHOSDigit * digit;
-    cout << "EMC sdigits " << endl ;
-    cout << "Digit 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->GetEntriesFast()) && 
-        (((AliPHOSDigit * )  sdigits->At(index))->GetId() <= maxEmc) ; index++) {
-      digit = (AliPHOSDigit * )  sdigits->At(index) ;
-      if(digit->GetNprimary() == 0) continue;
-      cout << setw(6)  <<  digit->GetId() << "   "  <<         setw(10)  <<  digit->GetAmp() <<   "    "  
-          << setw(6)  <<  digit->GetIndexInList() << "    "   
-          << setw(5)  <<  digit->GetNprimary() <<"    ";
-      
+        ((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) ) ; 
+      }  
     }    
-    cout << endl;
   }
 
   if(strstr(option,"all")||strstr(option,"CPV")){
     
     //loop over CPV digits
     AliPHOSDigit * digit;
-    cout << "CPV sdigits " << endl ;
-    cout << "Digit Id    Amplitude     Index     Nprim  Primaries list " <<  endl;      
+    printf("\nCPV sdigits\n") ; 
     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
     Int_t index ;
     for (index = 0 ; index < sdigits->GetEntriesFast(); index++) {
-      digit = (AliPHOSDigit * )  sdigits->At(index) ;
+      digit = dynamic_cast<AliPHOSDigit *>( sdigits->At(index) ) ;
       if(digit->GetId() > maxEmc){
-       cout << setw(6)  <<  digit->GetId() << "   "  <<        setw(10)  <<  digit->GetAmp() <<   "    "  
-            << setw(6)  <<  digit->GetIndexInList() << "    "   
-            << setw(5)  <<  digit->GetNprimary() <<"    ";
-       
+       printf("\n%6d  %8d    %4d      %2d :",
+               digit->GetId(), digit->GetAmp(), 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) ) ; 
+       }
       }    
     }
   }
-
 }
 
 //____________________________________________________________________________ 
-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() ; 
 }