AliRoot Core  ee782a0 (ee782a0)
MUONTriggerChamberEfficiency.C
Go to the documentation of this file.
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 #include "TMap.h"
28 #include "TSystem.h"
29 #include "TMath.h"
30 #include "TArrayI.h"
31 
32 // MUON includes
33 #include "AliMUONCDB.h"
34 #include "AliMUONCalibrationData.h"
37 #include "AliCDBManager.h"
38 #include "AliCDBRunRange.h"
39 #include "AliCDBId.h"
40 #include "AliCDBPath.h"
41 #include "AliCDBEntry.h"
42 
43 #endif
44 
53 
54 void MUONTriggerChamberEfficiency(TString inputFile = "./MUON.TriggerEfficiencyMap.root",
55  TString outputCDB = "",
56  Int_t firstRun=0, Int_t lastRun = AliCDBRunRange::Infinity()
57 )
58 {
68 
69  AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
70 
71  AliMUONTriggerEfficiencyCells* effMap = new AliMUONTriggerEfficiencyCells(inputFile.Data());
72 
73  if ( outputCDB.IsNull() ){
74  // Draw the efficiency and exit
75  AliCDBManager::Instance()->SetRun(firstRun);
77 
78  trigChEff->DisplayEfficiency(kFALSE,kFALSE);
79  return;
80  }
81 
82 
83  // Write efficiency on OCDB
84 
85  AliCDBManager::Instance()->SetSpecificStorage("MUON/Calib/TriggerEfficiency", outputCDB.Data());
86 
87  AliMUONCDB::WriteToCDB(effMap, "MUON/Calib/TriggerEfficiency", firstRun, lastRun, "Measured efficiencies");
88 }
89 
90 //____________________________________________________________
91 void ShowOCDBmap(Int_t runNumber = 0, TString specificCDB="", TString ocdbPath = "local://$ALICE_ROOT/OCDB", TString runType="Full")
92 {
99  if ( ocdbPath.BeginsWith("alien://") || ocdbPath.BeginsWith("raw://"))
100  TGrid::Connect("alien://");
101 
102  if (!ocdbPath.CompareTo("MC"))
103  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data(),runType.Data());
104  else
105  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath.Data());
106 
107  if ( !specificCDB.IsNull() )
108  AliCDBManager::Instance()->SetSpecificStorage("MUON/Calib/TriggerEfficiency", specificCDB.Data());
111 
113  trigChEff->DisplayEfficiency(kFALSE,kFALSE);
114 }
115 
116 
117 //____________________________________________________________
118 void FillHisto(TH1* histo, Int_t nevents)
119 {
121  for ( Int_t ibin=1; ibin<=histo->GetXaxis()->GetNbins(); ++ibin ) {
122  Double_t binCenter = histo->GetXaxis()->GetBinCenter(ibin);
123  for ( Int_t ievent=0; ievent<nevents; ++ievent ) {
124  histo->Fill(binCenter);
125  }
126  }
127 }
128 
129 //____________________________________________________________
130 void BuildDefaultMap(TString outFilename="/tmp/defTrigChEff.root", Double_t globalValue = 1., Int_t nevents = 100000)
131 {
133 
134  // Create histograms
135  enum { kBendingEff, kNonBendingEff, kBothPlanesEff, kAllTracks, kNcounts};
136  TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
137 
138  const Char_t* yAxisTitle = "counts";
139 
140  const Int_t kNboards = 234; //AliMpConstants::NofLocalBoards();
141  const Int_t kFirstTrigCh = 11;//AliMpConstants::NofTrackingChambers()+1;
142  const Int_t kNchambers = 4;
143  const Int_t kNslats = 18;
144 
145  Int_t chamberBins = kNchambers;
146  Float_t chamberLow = kFirstTrigCh-0.5, chamberHigh = kFirstTrigCh+kNchambers-0.5;
147  const Char_t* chamberName = "chamber";
148 
149  Int_t slatBins = kNslats;
150  Float_t slatLow = 0-0.5, slatHigh = kNslats-0.5;
151  const Char_t* slatName = "slat";
152 
153  Int_t boardBins = kNboards;
154  Float_t boardLow = 1-0.5, boardHigh = kNboards+1.-0.5;
155  const Char_t* boardName = "board";
156 
157  TString baseName, histoName, histoTitle;
158  TList* histoList = new TList();
159  histoList->SetOwner();
160 
161  TH1F* histo;
162 
163  for(Int_t icount=0; icount<kNcounts; icount++){
164  histoName = Form("%sCountChamber", countTypeName[icount].Data());
165  histo = new TH1F(histoName, histoName,
166  chamberBins, chamberLow, chamberHigh);
167  histo->GetXaxis()->SetTitle(chamberName);
168  histo->GetYaxis()->SetTitle(yAxisTitle);
169  Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
170  FillHisto(histo, (Int_t)nfills);
171  histoList->AddLast(histo);
172  } // loop on counts
173 
174  for(Int_t icount=0; icount<kNcounts; icount++){
175  for(Int_t ch=0; ch<kNchambers; ch++){
176  histoName = Form("%sCountSlatCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
177  histo = new TH1F(histoName, histoName,
178  slatBins, slatLow, slatHigh);
179  histo->GetXaxis()->SetTitle(slatName);
180  histo->GetYaxis()->SetTitle(yAxisTitle);
181  Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
182  FillHisto(histo, (Int_t)nfills);
183  histoList->AddLast(histo);
184  } // loop on chamber
185  } // loop on counts
186 
187  for(Int_t icount=0; icount<kNcounts; icount++){
188  for(Int_t ch=0; ch<kNchambers; ch++){
189  histoName = Form("%sCountBoardCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
190  histo = new TH1F(histoName, histoName,
191  boardBins, boardLow, boardHigh);
192  histo->GetXaxis()->SetTitle(boardName);
193  histo->GetYaxis()->SetTitle(yAxisTitle);
194  Double_t nfills = ( icount == kAllTracks ) ? nevents : globalValue * (Double_t)nevents;
195  FillHisto(histo, (Int_t)nfills);
196  histoList->AddLast(histo);
197  } // loop on chamber
198  } // loop on counts
199 
200  TFile* outFile = TFile::Open(outFilename,"create");
201  histoList->Write("triggerChamberEff",TObject::kSingleKey);
202  outFile->Close();
203 }
204 
205 
206 //____________________________________________________________
207 void BuildSystematicMap ( TString runList, Double_t globalSyst = -0.02, TString outputCDB = "local://CDB", TString inputCDB = "raw://", TString systematicPerBoard = "" )
208 {
214 
215  // Read run list
216  TArrayI runListArray(1000);
217  Int_t nruns=0;
218  if ( gSystem->AccessPathName(runList.Data()) ) {
219  runList.ReplaceAll(","," ");
220  TObjArray* arr = runList.Tokenize(" ");
221  TObjString* objStr = 0x0;
222  TIter nextRun(arr);
223  while ( (objStr = static_cast<TObjString*>(nextRun())) ) {
224  if ( objStr->String().IsDigit() ) runListArray[nruns++] = objStr->String().Atoi();
225  }
226  delete arr;
227  }
228  else {
229  ifstream inRunList(runList.Data());
230  TString currLine = "";
231  while ( ! inRunList.eof() ) {
232  currLine.ReadLine(inRunList);
233  currLine.ReplaceAll(","," ");
234  TObjArray* arr = currLine.Tokenize(" ");
235  TObjString* objStr = 0x0;
236  TIter nextRun(arr);
237  while ( (objStr = static_cast<TObjString*>(nextRun())) ) {
238  if ( objStr->String().IsDigit() ) runListArray[nruns++] = objStr->String().Atoi();
239  }
240  delete arr;
241  }
242  inRunList.close();
243  }
244  runListArray.Set(nruns);
245 
246  // Read maps from input CDB
248  mgr->SetDefaultStorage(inputCDB.Data());
249  TMap effMapsList;
250  AliCDBId* prevId = 0x0;
251  for ( Int_t irun=0; irun<runListArray.GetSize(); irun++ ) {
252  mgr->SetRun(runListArray[irun]);
253  AliCDBEntry* entry = mgr->Get(AliCDBPath("MUON/Calib/TriggerEfficiency"));
254  const AliCDBId cdbId = entry->GetId();
255  if ( prevId && cdbId.IsEqual(prevId) ) continue;
256  delete prevId;
257  prevId = static_cast<AliCDBId*>(cdbId.Clone());
258 
260  effMapsList.Add(new TObjString(Form("%i_%i",cdbId.GetFirstRun(),cdbId.GetLastRun())),effMap->Clone());
261  }
262  delete prevId;
263 
264 
265  // Set global systematic variation
266  Int_t nPoints = 234*4*3;
267  TArrayD systDiff(nPoints);
268  systDiff.Reset(globalSyst);
269  if ( ! systematicPerBoard.IsNull() && gSystem->AccessPathName(systematicPerBoard.Data()) == 0 ) {
270  // If file with specific chamber variation exists, read values
271  TString currLine = "";
272  ifstream inFile(systematicPerBoard.Data());
273  while ( ! inFile.eof() ) {
274  currLine.ReadLine(inFile);
275  TObjArray* arr = currLine.Tokenize(" ");
276  if ( arr->GetEntries() >= 4 && static_cast<TObjString*>(arr->UncheckedAt(0))->String().IsDigit() ) {
277  Int_t ich = static_cast<TObjString*>(arr->UncheckedAt(0))->String().Atoi()%11;
278  Int_t iboard = static_cast<TObjString*>(arr->UncheckedAt(1))->String().Atoi()-1;
279  Int_t icount = static_cast<TObjString*>(arr->UncheckedAt(2))->String().Atoi();
280  Double_t systErr = static_cast<TObjString*>(arr->UncheckedAt(3))->String().Atof();
281  Int_t idx = ich*234*3+iboard*3+icount;
282  if ( idx < nPoints ) systDiff[idx] = systErr;
283  }
284  delete arr;
285  }
286  inFile.close();
287 
288  // Set both planes uncertainty
289  // by default take the maximum uncertainty between bending and non-bending
290  for ( Int_t ich=0; ich<4; ich++ ) {
291  for ( Int_t iboard=0; iboard<234; iboard++ ) {
292  Int_t idx = ich*234*3+iboard*3;
293  Double_t diffBend = systDiff[idx];
294  Double_t diffNonBend = systDiff[idx+1];
295  systDiff[idx+2] = ( TMath::Abs(diffBend) > TMath::Abs(diffNonBend) ) ? diffBend : diffNonBend;
296  }
297  }
298  }
299 
300  mgr->SetDefaultStorage(outputCDB.Data());
301 
302  enum { kBendingEff, kNonBendingEff, kBothPlanesEff};
303  TString countTypeName[3] = {"bendPlane", "nonBendPlane","bothPlanes"};
304 
305  TIter next(&effMapsList);
306  TObject* key = 0x0;
307  while ( ( key = next() ) ) {
308  TString runRange = key->GetName();
309  AliMUONTriggerEfficiencyCells* effMap = static_cast<AliMUONTriggerEfficiencyCells*>(effMapsList.GetValue(key));
310  for ( Int_t ich=0; ich<4; ich++ ) {
311  TString histoName = Form("allTracksCountBoardCh%i", 11+ich);
312  TH1* auxHisto = static_cast<TH1*>(effMap->GetHistoList()->FindObject(histoName));
313  for ( Int_t icount=0; icount<3; icount++ ) {
314  histoName = Form("%sCountBoardCh%i", countTypeName[icount].Data(), 11+ich);
315  TH1* histo = static_cast<TH1*>(effMap->GetHistoList()->FindObject(histoName));
316  for ( Int_t iboard=0; iboard<234; iboard++ ) {
317  Int_t idx = ich*234*3+iboard*3+icount;
318  Int_t ibin = iboard+1;
319  Double_t countAll = auxHisto->GetBinContent(ibin);
320  Double_t countDiff = TMath::Nint(systDiff[idx] * countAll);
321  Double_t countOrig = histo->GetBinContent(ibin);
322  Double_t newCount = countOrig+countDiff;
323  if ( newCount < 0 || newCount > countAll ) {
324  printf("WARNING: ch %i board %i cath %i systDiff %g newEff %g / %g\n",11+ich,ibin,icount,systDiff[idx],newCount,countAll);
325  if ( newCount < 0 ) newCount = 0;
326  else newCount = countAll;
327  printf(" => setting numerator to %g\n",newCount);
328  }
329  histo->SetBinContent(ibin,newCount);
330  if ( histo->GetSumw2N() > 0 ) histo->SetBinError(ibin,TMath::Sqrt(newCount));
331  }
332  }
333  }
334  TObjArray* arr = runRange.Tokenize("_");
335  Int_t firstRun = static_cast<TObjString*>(arr->UncheckedAt(0))->String().Atoi();
336  Int_t lastRun = static_cast<TObjString*>(arr->UncheckedAt(1))->String().Atoi();
337  AliMUONCDB::WriteToCDB(effMap, "MUON/Calib/TriggerEfficiency", firstRun, lastRun, "Measured efficiencies");
338  }
339 }
340 
341 //____________________________________________________________
342 void CompleteEfficiency(TString effFileWithHoles, TString effFileCompatible, TString outFilename)
343 {
344  //
347  //
348 
349  TList* histoList[2] = {0x0, 0x0};
350  TString filenames[2] = {effFileWithHoles, effFileCompatible};
351  for ( Int_t ifile=0; ifile<2; ifile++ ) {
352  TFile* file = TFile::Open(filenames[ifile].Data());
353  if ( ! file ) {
354  printf("Fatal: cannot find %s\n", filenames[ifile].Data());
355  return;
356  }
357  histoList[ifile] = static_cast<TList*> (file->FindObjectAny("triggerChamberEff"));
358  if ( ! histoList[ifile] ) {
359  printf("Cannot find histo list in file %s\n", filenames[ifile].Data());
360  return;
361  }
362  }
363 
364  TString detElemName[2] = {"Slat", "Board"};
365  enum { kBendingEff, kNonBendingEff, kBothPlanesEff, kAllTracks, kNcounts};
366  TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
367 
368  Bool_t isChanged = kFALSE;
369  TString histoName = "";
370  for ( Int_t idet=0; idet<2; idet++ ) {
371  for ( Int_t ich=11; ich<=14; ich++ ) {
372  histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), ich);
373  TH1* allTracksHisto = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
374  for ( Int_t ibin=1; ibin<=allTracksHisto->GetXaxis()->GetNbins(); ibin++ ) {
375  if ( allTracksHisto->GetBinContent(ibin) > 0. ) continue;
376  isChanged = kTRUE;
377  printf("Modifying info for Ch %i %s %3i\n", ich, detElemName[idet].Data(), (Int_t)allTracksHisto->GetXaxis()->GetBinCenter(ibin));
378  // If allTracks has no entries, it means that efficiency could not be calculated for this bin:
379  // fill information from the compatible histogram
380 
381  // Check the statistics collected by the switched off detection element
382  Double_t nTracks = 0;
383  for ( Int_t jch=11; jch<=14; jch++ ) {
384  histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), jch);
385  TH1* allTracksOtherCh = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
386  nTracks = allTracksOtherCh->GetBinContent(ibin);
387  if ( nTracks > 0. ) {
388  //printf("Statistics for %s : %g\n", histoName.Data(), nTracks); // REMEMBER TO CUT
389  break;
390  }
391  }
392 
393  histoName = Form("%sCount%sCh%i", countTypeName[kAllTracks].Data(), detElemName[idet].Data(), ich);
394  TH1* allTracksHistoAux = static_cast<TH1*> (histoList[1]->FindObject(histoName.Data()));
395  Double_t nTracksNew = allTracksHistoAux->GetBinContent(ibin);
396  if ( nTracksNew == 0.) {
397  printf("Warning: new histogram has no entries for Ch %i %s %3i\n", ich, detElemName[idet].Data(), (Int_t)allTracksHisto->GetXaxis()->GetBinCenter(ibin));
398  continue;
399  }
400  Double_t scaleFactor = TMath::Min(nTracksNew, nTracks) / nTracksNew;
401  //printf("Statistics ineff %g new %g scaleFactor %g\n", nTracks, nTracksNew, scaleFactor); // REMEMBER TO CUT
402 
403  for ( Int_t icount=0; icount<kNcounts; icount++ ) {
404  histoName = Form("%sCount%sCh%i", countTypeName[icount].Data(), detElemName[idet].Data(), ich);
405  TH1* auxHisto = static_cast<TH1*> (histoList[1]->FindObject(histoName.Data()));
406  TH1* histo = static_cast<TH1*> (histoList[0]->FindObject(histoName.Data()));
407  histo->SetBinContent(ibin, auxHisto->GetBinContent(ibin) * scaleFactor);
408  if ( histo->GetSumw2N() > 0 ) histo->SetBinError(ibin, auxHisto->GetBinError(ibin)*scaleFactor);
409  } // loop on cont types
410  } // loop on histogram bins
411  } // loop on chamber
412  } // loop on detection element (slat or board)
413 
414  if ( ! isChanged ) {
415  printf("Input histograms not modified\n");
416  return;
417  }
418  TFile* outFile = TFile::Open(outFilename,"create");
419  histoList[0]->Write("triggerChamberEff",TObject::kSingleKey);
420  outFile->Close();
421 }
422 
423 
void ShowOCDBmap(Int_t runNumber=0, TString specificCDB="", TString ocdbPath="local://$ALICE_ROOT/OCDB", TString runType="Full")
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TFile * Open(const char *filename, Long64_t &nevents)
#define TObjArray
void FillHisto(TH1 *histo, Int_t nevents)
TList * GetHistoList()
Get list of histograms.
Store and give access to the trigger chamber efficiency.
Int_t GetLastRun() const
Definition: AliCDBId.h:49
AliCDBEntry * Get(const AliCDBId &query, Bool_t forceCaching=kFALSE)
TObject * GetObject()
Definition: AliCDBEntry.h:56
void SetSpecificStorage(const char *calibType, const char *dbString, Int_t version=-1, Int_t subVersion=-1)
virtual Bool_t IsEqual(const TObject *obj) const
Definition: AliCDBId.cxx:147
void MUONTriggerChamberEfficiency(TString inputFile="./MUON.TriggerEfficiencyMap.root", TString outputCDB="", Int_t firstRun=0, Int_t lastRun=AliCDBRunRange::Infinity())
void CompleteEfficiency(TString effFileWithHoles, TString effFileCompatible, TString outFilename)
void SetRun(Int_t run)
Definition: AliCDBEntry.h:18
AliMUONTriggerEfficiencyCells * TriggerEfficiency() const
Get the trigger efficiency map.
void SetDefaultStorage(const char *dbString)
void DisplayEfficiency(Bool_t perSlat=kFALSE, Bool_t show2Dhisto=kTRUE)
void WriteToCDB(const char *calibpath, TObject *object, Int_t startRun, Int_t endRun, const char *filename)
Single entry point to access MUON calibration data.
void BuildDefaultMap(TString outFilename="/tmp/defTrigChEff.root", Double_t globalValue=1., Int_t nevents=100000)
static Int_t runNumber
Definition: pdc06_config.C:126
Calculate, apply and possibly draw trigger chamber efficiency.
AliCDBId & GetId()
Definition: AliCDBEntry.h:51
static AliCDBManager * Instance(TMap *entryCache=NULL, Int_t run=-1)
void BuildSystematicMap(TString runList, Double_t globalSyst=-0.02, TString outputCDB="local://CDB", TString inputCDB="raw://", TString systematicPerBoard="")
Int_t GetFirstRun() const
Definition: AliCDBId.h:48
static Int_t Infinity()