MUON/macros/MUONTriggerChamberEfficiency.C

   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 #if !defined(__CINT__) || defined(__MAKECINT__)
  19
  20 #include "Riostream.h"
  21
  22 // ROOT includes
  23 #include "TGrid.h"
  24 #include "TString.h"
  25 #include "TFile.h"
  26 #include "TH1.h"
  27
  28 // MUON includes
  29 #include "AliMUONCDB.h"
  30 #include "AliMUONCalibrationData.h"
  31 #include "AliMUONTriggerEfficiencyCells.h"
  32 #include "AliMUONTriggerChamberEfficiency.h"
  33 #include "AliCDBManager.h"
  34 #include "AliCDBRunRange.h"
  35
  36 #endif
  37
  38 /// \ingroup macros
  39 /// \file MUONTriggerChamberEfficiency.C
  40 /// \brief Macro to view and save the trigger chamber efficiency map
  41 /// calculated during reconstruction.
  42 ///
  43 /// Efficiency map can be made available for next simulation.
  44 ///
  45 /// \author Diego Stocco, Subatech, Nantes
  46
  47 void MUONTriggerChamberEfficiency(TString inputFile = "./MUON.TriggerEfficiencyMap.root",
  48                                   TString outputCDB = "",
  49                                   Int_t firstRun=0, Int_t lastRun = AliCDBRunRange::Infinity()
  50 )
  51 {
  52 /// \param inputFile (default "./MUON.TriggerEfficiencyMaps.root")
  53 ///     File with the numerator and denominator histos for efficiency calculation
  54 ///     (It is the output of the PWG3/muon/AliAnalysisTaskTrigChEff analysis
  55 /// \param outputCDB (default "")
  56 ///     add the map on the specified CDB
  57 /// \param firstRun (default 0)
  58 ///     first run of validity for CDB object
  59 /// \param lastRun (default AliCDBRunRange::Infinity())
  60 ///     last run of validity for CDB Object
  61
  62   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  63
  64   AliMUONTriggerEfficiencyCells* effMap = new AliMUONTriggerEfficiencyCells(inputFile.Data());
  65
  66   if ( outputCDB.IsNull() ){
  67     // Draw the efficiency and exit
  68     AliCDBManager::Instance()->SetRun(firstRun);
  69     AliMUONTriggerChamberEfficiency* trigChEff = new AliMUONTriggerChamberEfficiency(effMap);
  70
  71     trigChEff->DisplayEfficiency(kFALSE,kFALSE);
  72     return;
  73   }
  74
  75
  76   // Write efficiency on OCDB
  77
  78   AliCDBManager::Instance()->SetSpecificStorage("MUON/Calib/TriggerEfficiency", outputCDB.Data());
  79
  80   AliMUONCDB::WriteToCDB(effMap, "MUON/Calib/TriggerEfficiency", firstRun, lastRun, "Measured efficiencies");
  81 }
  82
  83 //____________________________________________________________
  84 void ShowOCDBmap(Int_t runNumber = 0, TString specificCDB="", TString ocdbPath = "local://$ALICE_ROOT/OCDB", TString runType="Full")
  85 {
  86 /// \param runNumber (default 0)
  87 ///     run number
  88 /// \param specificCDB (default "")
  89 ///     specific CDB for trigger efficiency
  90 /// \param ocdbPath(default "local://$ALICE_ROOT/OCDB")
  91 ///     path to OCDB
  92   if ( ocdbPath.BeginsWith("alien://") || ocdbPath.BeginsWith("raw://"))
  93     TGrid::Connect("alien://");
  94
  95   if (!ocdbPath.CompareTo("MC"))
  96     AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data(),runType.Data());
  97   else
  98     AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data());
  99
 100   if ( !specificCDB.IsNull() )
 101     AliCDBManager::Instance()->SetSpecificStorage("MUON/Calib/TriggerEfficiency", specificCDB.Data());
 102   AliCDBManager::Instance()->SetRun(runNumber);
 103   AliMUONCalibrationData calib(runNumber);
 104
 105   AliMUONTriggerChamberEfficiency* trigChEff = new AliMUONTriggerChamberEfficiency(calib.TriggerEfficiency());
 106   trigChEff->DisplayEfficiency(kFALSE,kFALSE);
 107 }
 108
 109
 110 //____________________________________________________________
 111 void FillHisto(TH1* histo, Int_t nevents)
 112 {
 113   /// Fill histogram with global value
 114   for ( Int_t ibin=1; ibin<=histo->GetXaxis()->GetNbins(); ++ibin ) {
 115     Double_t binCenter = histo->GetXaxis()->GetBinCenter(ibin);
 116     for ( Int_t ievent=0; ievent<nevents; ++ievent ) {
 117       histo->Fill(binCenter);
 118     }
 119   }
 120 }
 121
 122 //____________________________________________________________
 123 void BuildDefaultMap(TString outFilename="/tmp/defTrigChEff.root", Double_t globalValue = 1., Int_t nevents = 100000)
 124 {
 125   /// Build default map (all boards with the same chosen value)
 126
 127   // Create histograms
 128   enum { kBendingEff, kNonBendingEff, kBothPlanesEff, kAllTracks, kNcounts};
 129   TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
 130
 131   const Char_t* yAxisTitle = "counts";
 132
 133   const Int_t kNboards = 234; //AliMpConstants::NofLocalBoards();
 134   const Int_t kFirstTrigCh = 11;//AliMpConstants::NofTrackingChambers()+1;
 135   const Int_t kNchambers = 4;
 136   const Int_t kNslats = 18;
 137
 138   Int_t chamberBins = kNchambers;
 139   Float_t chamberLow = kFirstTrigCh-0.5, chamberHigh = kFirstTrigCh+kNchambers-0.5;
 140   const Char_t* chamberName = "chamber";
 141
 142   Int_t slatBins = kNslats;
 143   Float_t slatLow = 0-0.5, slatHigh = kNslats-0.5;
 144   const Char_t* slatName = "slat";
 145
 146   Int_t boardBins = kNboards;
 147   Float_t boardLow = 1-0.5, boardHigh = kNboards+1.-0.5;
 148   const Char_t* boardName = "board";
 149
 150   TString baseName, histoName, histoTitle;
 151   TList* histoList = new TList();
 152   histoList->SetOwner();
 153
 154   TH1F* histo;
 155
 156   for(Int_t icount=0; icount<kNcounts; icount++){
 157     histoName = Form("%sCountChamber", countTypeName[icount].Data());
 158     histo = new TH1F(histoName, histoName,
 159                      chamberBins, chamberLow, chamberHigh);
 160     histo->GetXaxis()->SetTitle(chamberName);
 161     histo->GetYaxis()->SetTitle(yAxisTitle);
 162     Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
 163     FillHisto(histo, (Int_t)nfills);
 164     histoList->AddLast(histo);
 165   } // loop on counts
 166
 167   for(Int_t icount=0; icount<kNcounts; icount++){
 168     for(Int_t ch=0; ch<kNchambers; ch++){
 169       histoName = Form("%sCountSlatCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
 170       histo = new TH1F(histoName, histoName,
 171                        slatBins, slatLow, slatHigh);
 172       histo->GetXaxis()->SetTitle(slatName);
 173       histo->GetYaxis()->SetTitle(yAxisTitle);
 174       Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
 175       FillHisto(histo, (Int_t)nfills);
 176       histoList->AddLast(histo);
 177     } // loop on chamber
 178   } // loop on counts
 179
 180   for(Int_t icount=0; icount<kNcounts; icount++){
 181     for(Int_t ch=0; ch<kNchambers; ch++){
 182       histoName = Form("%sCountBoardCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
 183       histo = new TH1F(histoName, histoName,
 184                        boardBins, boardLow, boardHigh);
 185       histo->GetXaxis()->SetTitle(boardName);
 186       histo->GetYaxis()->SetTitle(yAxisTitle);
 187       Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
 188       FillHisto(histo, (Int_t)nfills);
 189       histoList->AddLast(histo);
 190     } // loop on chamber
 191   } // loop on counts
 192
 193   TFile* outFile = TFile::Open(outFilename,"create");
 194   histoList->Write("triggerChamberEff",TObject::kSingleKey);
 195   outFile->Close();
 196 }
 197
 198
 199 //____________________________________________________________
 200 void CompleteEfficiency(TString effFileWithHoles, TString effFileCompatible, TString outFilename)
 201 {
 202   //
 203   /// When a local board or RPC is missing, the efficiency of other boards cannot be calculated
 204   /// If an efficiency file of the same period is available, it could be used to fill the missing information
 205   //
 206
 207   TList* histoList[2] = {0x0, 0x0};
 208   TString filenames[2] = {effFileWithHoles, effFileCompatible};
 209   for ( Int_t ifile=0; ifile<2; ifile++ ) {
 210     TFile* file = TFile::Open(filenames[ifile].Data());
 211     if ( ! file ) {
 212       printf("Fatal: cannot find %s\n", filenames[ifile].Data());
 213       return;
 214     }
 215     histoList[ifile] = static_cast<TList*> (file->FindObjectAny("triggerChamberEff"));
 216     if ( ! histoList[ifile] ) {
 217       printf("Cannot find histo list in file %s\n", filenames[ifile].Data());
 218       return;
 219     }
 220   }
 221
 222   TString detElemName[2] = {"Slat", "Board"};
 223   enum { kBendingEff, kNonBendingEff, kBothPlanesEff, kAllTracks, kNcounts};
 224   TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
 225
 226   Bool_t isChanged = kFALSE;
 227   TString histoName = "";
 228   for ( Int_t idet=0; idet<2; idet++ ) {
 229     for ( Int_t ich=11; ich<=14; ich++ ) {
 230       histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), ich);
 231       TH1* allTracksHisto = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
 232       for ( Int_t ibin=1; ibin<=allTracksHisto->GetXaxis()->GetNbins(); ibin++ ) {
 233         if ( allTracksHisto->GetBinContent(ibin) > 0. ) continue;
 234         isChanged = kTRUE;
 235         printf("Modifying info for Ch %i %s %3i\n", ich, detElemName[idet].Data(), (Int_t)allTracksHisto->GetXaxis()->GetBinCenter(ibin));
 236         // If allTracks has no entries, it means that efficiency could not be calculated for this bin:
 237         // fill information from the compatible histogram
 238
 239         // Check the statistics collected by the switched off detection element
 240         Double_t nTracks = 0;
 241         for ( Int_t jch=11; jch<=14; jch++ ) {
 242           histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), jch);
 243           TH1* allTracksOtherCh = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
 244           nTracks = allTracksOtherCh->GetBinContent(ibin);
 245           if ( nTracks > 0. ) {
 246             //printf("Statistics for %s : %g\n", histoName.Data(), nTracks); // REMEMBER TO CUT
 247             break;
 248           }
 249         }
 250
 251         histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), ich);
 252         TH1* allTracksHistoAux = static_cast<TH1*> (histoList[1]->FindObject(histoName.Data()));
 253         Double_t nTracksNew = allTracksHistoAux->GetBinContent(ibin);
 254         if ( nTracksNew == 0.) {
 255           printf("Warning: new histogram has no entries for Ch %i %s %3i\n", ich, detElemName[idet].Data(), (Int_t)allTracksHisto->GetXaxis()->GetBinCenter(ibin));
 256           continue;
 257         }
 258         Double_t scaleFactor = TMath::Min(nTracksNew, nTracks) / nTracksNew;
 259         //printf("Statistics ineff %g  new %g  scaleFactor %g\n", nTracks, nTracksNew, scaleFactor); // REMEMBER TO CUT
 260
 261         for ( Int_t icount=0; icount<kNcounts; icount++ ) {
 262           histoName = Form("%sCount%sCh%i", countTypeName[icount].Data(), detElemName[idet].Data(), ich);
 263           TH1* auxHisto = static_cast<TH1*> (histoList[1]->FindObject(histoName.Data()));
 264           TH1* histo = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
 265           histo->SetBinContent(ibin, auxHisto->GetBinContent(ibin) * scaleFactor);
 266           if ( histo->GetSumw2N() > 0 ) histo->SetBinError(ibin, auxHisto->GetBinError(ibin)*scaleFactor);
 267         } // loop on cont types
 268       } // loop on histogram bins
 269     } // loop on chamber
 270   } // loop on detection element (slat or board)
 271
 272   if ( ! isChanged ) {
 273     printf("Input histograms not modified\n");
 274     return;
 275   }
 276   TFile* outFile = TFile::Open(outFilename,"create");
 277   histoList[0]->Write("triggerChamberEff",TObject::kSingleKey);
 278   outFile->Close();
 279 }
 280
 281