Flexible pt range for the efficiency histogramming
[u/mrichter/AliRoot.git] / MUON / 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