Fixed memory leak
[u/mrichter/AliRoot.git] / MUON / MUONTriggerChamberEfficiency.C
index 95bda5c..61a7d99 100644 (file)
 /* $Id$ */
 
 #if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "Riostream.h"
+
 // ROOT includes
 #include "TGrid.h"
 #include "TString.h"
+#include "TFile.h"
+#include "TH1.h"
 
 // MUON includes
 #include "AliMUONCDB.h"
@@ -27,7 +32,6 @@
 #include "AliMUONTriggerChamberEfficiency.h"
 #include "AliCDBManager.h"
 #include "AliCDBRunRange.h"
-#include "Riostream.h"
 
 #endif
 
 ///
 /// \author Diego Stocco, Subatech, Nantes
 
-void MUONTriggerChamberEfficiency(const Char_t* inputFile="./MUON.TriggerEfficiencyMap.root",
-                                 Bool_t addMapInLocalOCDB = kFALSE,
-                                 Int_t firstRun=0, Int_t lastRun = AliCDBRunRange::Infinity())
+void MUONTriggerChamberEfficiency(TString inputFile = "./MUON.TriggerEfficiencyMap.root",
+                                 TString outputCDB = "",
+                                 Int_t firstRun=0, Int_t lastRun = AliCDBRunRange::Infinity()
+)
 {
 /// \param inputFile (default "./MUON.TriggerEfficiencyMaps.root")
 ///     File with the numerator and denominator histos for efficiency calculation
 ///     (It is the output of the PWG3/muon/AliAnalysisTaskTrigChEff analysis
-/// \param addMapInLocalOCDB (default kFALSE)
-///     if the map is OK, add it in local OCDB
+/// \param outputCDB (default "")
+///     add the map on the specified CDB
 /// \param firstRun (default 0)
 ///     first run of validity for CDB object
 /// \param lastRun (default AliCDBRunRange::Infinity())
@@ -56,32 +61,221 @@ void MUONTriggerChamberEfficiency(const Char_t* inputFile="./MUON.TriggerEfficie
 
   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 
-  AliMUONTriggerEfficiencyCells* effMap = new AliMUONTriggerEfficiencyCells(inputFile);
-  if ( ! addMapInLocalOCDB ){
+  AliMUONTriggerEfficiencyCells* effMap = new AliMUONTriggerEfficiencyCells(inputFile.Data());
+
+  if ( outputCDB.IsNull() ){
+    // Draw the efficiency and exit
     AliCDBManager::Instance()->SetRun(firstRun);
-    AliMUONTriggerChamberEfficiency trigChEff(effMap);
+    AliMUONTriggerChamberEfficiency* trigChEff = new AliMUONTriggerChamberEfficiency(effMap);
  
-    trigChEff.DisplayEfficiency();
+    trigChEff->DisplayEfficiency(kFALSE,kFALSE);
+    return;
   }
-  else
-    AliMUONCDB::WriteToCDB("MUON/Calib/TriggerEfficiency", effMap, firstRun, lastRun,true);
+
+
+  // Write efficiency on OCDB
+
+  AliCDBManager::Instance()->SetSpecificStorage("MUON/Calib/TriggerEfficiency", outputCDB.Data());
+  
+  AliMUONCDB::WriteToCDB(effMap, "MUON/Calib/TriggerEfficiency", firstRun, lastRun, "Measured efficiencies");
 }
 
 //____________________________________________________________
-void ShowOCDBmap(Int_t runNumber = 0, TString ocdbPath="local://$ALICE_ROOT/OCDB")
+void ShowOCDBmap(Int_t runNumber = 0, TString specificCDB="", TString ocdbPath = "local://$ALICE_ROOT/OCDB", TString runType="Full")
 {
 /// \param runNumber (default 0)
 ///     run number
+/// \param specificCDB (default "")
+///     specific CDB for trigger efficiency
 /// \param ocdbPath(default "local://$ALICE_ROOT/OCDB")
 ///     path to OCDB
   if ( ocdbPath.BeginsWith("alien://") || ocdbPath.BeginsWith("raw://"))
     TGrid::Connect("alien://");
 
-  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data());
+  if (!ocdbPath.CompareTo("MC"))
+    AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data(),runType.Data());
+  else
+    AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data());
+  
+  if ( !specificCDB.IsNull() )
+    AliCDBManager::Instance()->SetSpecificStorage("MUON/Calib/TriggerEfficiency", specificCDB.Data());
   AliCDBManager::Instance()->SetRun(runNumber);
   AliMUONCalibrationData calib(runNumber);
 
-  AliMUONTriggerChamberEfficiency trigChEff(calib.TriggerEfficiency());
-  trigChEff.DisplayEfficiency();
+  AliMUONTriggerChamberEfficiency* trigChEff = new AliMUONTriggerChamberEfficiency(calib.TriggerEfficiency());
+  trigChEff->DisplayEfficiency(kFALSE,kFALSE);
+}
+
+
+//____________________________________________________________
+void FillHisto(TH1* histo, Int_t nevents)
+{
+  /// Fill histogram with global value
+  for ( Int_t ibin=1; ibin<=histo->GetXaxis()->GetNbins(); ++ibin ) {
+    Double_t binCenter = histo->GetXaxis()->GetBinCenter(ibin);
+    for ( Int_t ievent=0; ievent<nevents; ++ievent ) {
+      histo->Fill(binCenter);
+    }
+  }
+}
+
+//____________________________________________________________
+void BuildDefaultMap(TString outFilename="/tmp/defTrigChEff.root", Double_t globalValue = 1., Int_t nevents = 100000)
+{
+  /// Build default map (all boards with the same chosen value)
+  
+  // Create histograms
+  enum { kBendingEff, kNonBendingEff, kBothPlanesEff, kAllTracks, kNcounts};
+  TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
+
+  const Char_t* yAxisTitle = "counts";
+
+  const Int_t kNboards = 234; //AliMpConstants::NofLocalBoards();
+  const Int_t kFirstTrigCh = 11;//AliMpConstants::NofTrackingChambers()+1;
+  const Int_t kNchambers = 4;
+  const Int_t kNslats = 18;
+
+  Int_t chamberBins = kNchambers;
+  Float_t chamberLow = kFirstTrigCh-0.5, chamberHigh = kFirstTrigCh+kNchambers-0.5;
+  const Char_t* chamberName = "chamber";
+
+  Int_t slatBins = kNslats;
+  Float_t slatLow = 0-0.5, slatHigh = kNslats-0.5;
+  const Char_t* slatName = "slat";
+
+  Int_t boardBins = kNboards;
+  Float_t boardLow = 1-0.5, boardHigh = kNboards+1.-0.5;
+  const Char_t* boardName = "board";
+
+  TString baseName, histoName, histoTitle;
+  TList* histoList = new TList();
+  histoList->SetOwner();
+  
+  TH1F* histo;
+
+  for(Int_t icount=0; icount<kNcounts; icount++){
+    histoName = Form("%sCountChamber", countTypeName[icount].Data());
+    histo = new TH1F(histoName, histoName,
+                     chamberBins, chamberLow, chamberHigh);
+    histo->GetXaxis()->SetTitle(chamberName);
+    histo->GetYaxis()->SetTitle(yAxisTitle);
+    Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
+    FillHisto(histo, (Int_t)nfills);
+    histoList->AddLast(histo);
+  } // loop on counts
+
+  for(Int_t icount=0; icount<kNcounts; icount++){
+    for(Int_t ch=0; ch<kNchambers; ch++){
+      histoName = Form("%sCountSlatCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
+      histo = new TH1F(histoName, histoName,
+                       slatBins, slatLow, slatHigh);
+      histo->GetXaxis()->SetTitle(slatName);
+      histo->GetYaxis()->SetTitle(yAxisTitle);
+      Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
+      FillHisto(histo, (Int_t)nfills);
+      histoList->AddLast(histo);
+    } // loop on chamber
+  } // loop on counts
+
+  for(Int_t icount=0; icount<kNcounts; icount++){
+    for(Int_t ch=0; ch<kNchambers; ch++){
+      histoName = Form("%sCountBoardCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
+      histo = new TH1F(histoName, histoName,
+                       boardBins, boardLow, boardHigh);
+      histo->GetXaxis()->SetTitle(boardName);
+      histo->GetYaxis()->SetTitle(yAxisTitle);
+      Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
+      FillHisto(histo, (Int_t)nfills);
+      histoList->AddLast(histo);
+    } // loop on chamber
+  } // loop on counts
+
+  TFile* outFile = TFile::Open(outFilename,"create");
+  histoList->Write("triggerChamberEff",TObject::kSingleKey);
+  outFile->Close();
 }
 
+
+//____________________________________________________________
+void CompleteEfficiency(TString effFileWithHoles, TString effFileCompatible, TString outFilename)
+{
+  //
+  /// When a local board or RPC is missing, the efficiency of other boards cannot be calculated
+  /// If an efficiency file of the same period is available, it could be used to fill the missing information
+  //
+  
+  TList* histoList[2] = {0x0, 0x0};
+  TString filenames[2] = {effFileWithHoles, effFileCompatible};
+  for ( Int_t ifile=0; ifile<2; ifile++ ) {
+    TFile* file = TFile::Open(filenames[ifile].Data());
+    if ( ! file ) {
+      printf("Fatal: cannot find %s\n", filenames[ifile].Data());
+      return;
+    }
+    histoList[ifile] = static_cast<TList*> (file->FindObjectAny("triggerChamberEff"));
+    if ( ! histoList[ifile] ) {
+      printf("Cannot find histo list in file %s\n", filenames[ifile].Data());
+      return;
+    }
+  }
+  
+  TString detElemName[2] = {"Slat", "Board"};
+  enum { kBendingEff, kNonBendingEff, kBothPlanesEff, kAllTracks, kNcounts};
+  TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
+  
+  Bool_t isChanged = kFALSE;
+  TString histoName = "";
+  for ( Int_t idet=0; idet<2; idet++ ) {
+    for ( Int_t ich=11; ich<=14; ich++ ) {
+      histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), ich);
+      TH1* allTracksHisto = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
+      for ( Int_t ibin=1; ibin<=allTracksHisto->GetXaxis()->GetNbins(); ibin++ ) {
+        if ( allTracksHisto->GetBinContent(ibin) > 0. ) continue;
+        isChanged = kTRUE;
+        printf("Modifying info for Ch %i %s %3i\n", ich, detElemName[idet].Data(), (Int_t)allTracksHisto->GetXaxis()->GetBinCenter(ibin));
+        // If allTracks has no entries, it means that efficiency could not be calculated for this bin:
+        // fill information from the compatible histogram
+        
+        // Check the statistics collected by the switched off detection element
+        Double_t nTracks = 0;
+        for ( Int_t jch=11; jch<=14; jch++ ) {
+          histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), jch);
+          TH1* allTracksOtherCh = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
+          nTracks = allTracksOtherCh->GetBinContent(ibin);
+          if ( nTracks > 0. ) {
+            //printf("Statistics for %s : %g\n", histoName.Data(), nTracks); // REMEMBER TO CUT
+            break;
+          }
+        }
+        
+        histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), ich);
+        TH1* allTracksHistoAux = static_cast<TH1*> (histoList[1]->FindObject(histoName.Data()));
+        Double_t nTracksNew = allTracksHistoAux->GetBinContent(ibin);
+        if ( nTracksNew == 0.) {
+          printf("Warning: new histogram has no entries for Ch %i %s %3i\n", ich, detElemName[idet].Data(), (Int_t)allTracksHisto->GetXaxis()->GetBinCenter(ibin));
+          continue;
+        }
+        Double_t scaleFactor = TMath::Min(nTracksNew, nTracks) / nTracksNew;
+        //printf("Statistics ineff %g  new %g  scaleFactor %g\n", nTracks, nTracksNew, scaleFactor); // REMEMBER TO CUT
+        
+        for ( Int_t icount=0; icount<kNcounts; icount++ ) {
+          histoName = Form("%sCount%sCh%i", countTypeName[icount].Data(), detElemName[idet].Data(), ich);
+          TH1* auxHisto = static_cast<TH1*> (histoList[1]->FindObject(histoName.Data()));
+          TH1* histo = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
+          histo->SetBinContent(ibin, auxHisto->GetBinContent(ibin) * scaleFactor);
+          if ( histo->GetSumw2N() > 0 ) histo->SetBinError(ibin, auxHisto->GetBinError(ibin)*scaleFactor);
+        } // loop on cont types
+      } // loop on histogram bins
+    } // loop on chamber
+  } // loop on detection element (slat or board)
+  
+  if ( ! isChanged ) {
+    printf("Input histograms not modified\n");
+    return;
+  }
+  TFile* outFile = TFile::Open(outFilename,"create");
+  histoList[0]->Write("triggerChamberEff",TObject::kSingleKey);
+  outFile->Close();  
+}
+
+