MUON + CheckCompiler
[u/mrichter/AliRoot.git] / MUON / MUONcalib / AliMUONTriggerChamberEfficiency.cxx
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 #include "AliMUONTriggerEfficiencyCells.h"
19 #include "AliMpConstants.h"
20
21 // Classes for display
22 #include "AliMUONTriggerDisplay.h"
23 #include "AliCDBManager.h"
24 #include "AliMpDDLStore.h"
25
26 #include "AliLog.h"
27
28 #include "TRandom.h"
29 #include "Riostream.h"
30 #include "TH1F.h"
31 #include "TObjArray.h"
32 #include "TGraphAsymmErrors.h"
33
34 #include "TH2F.h"
35 #include "TCanvas.h"
36 #include "TROOT.h"
37
38 #include "AliMUONTriggerChamberEfficiency.h"
39
40 //-----------------------------------------------------------------------------
41 /// \class AliMUONTriggerChamberEfficiency
42 /// A class to store and give access to the trigger chamber efficiency.
43 ///
44 /// Efficiency is stored per cathode on local boards
45 ///
46 /// The main method of this class is IsTriggered().
47 ///
48 /// \author Diego Stocco; INFN Torino
49 //-----------------------------------------------------------------------------
50
51 /// \cond CLASSIMP
52 ClassImp(AliMUONTriggerChamberEfficiency)
53 /// \endcond
54
55 //__________________________________________________________________________
56 AliMUONTriggerChamberEfficiency::AliMUONTriggerChamberEfficiency(AliMUONTriggerEfficiencyCells* effCells)
57 :
58 TObject(),
59 fIsOwner(kFALSE),
60 fEfficiencyMap(effCells),
61 fEfficiencyObjects(0x0),
62 fDisplayList(0x0)
63 {
64 ///  Default constructor.
65   FillFromList();
66 }
67
68 //__________________________________________________________________________
69 AliMUONTriggerChamberEfficiency::AliMUONTriggerChamberEfficiency(const Char_t* filename, const Char_t* listname)
70 :
71 TObject(),
72 fIsOwner(kTRUE),
73 fEfficiencyMap(0x0),
74 fEfficiencyObjects(0x0),
75 fDisplayList(0x0)
76 {
77 ///  Constructor using an ASCII file.
78   fEfficiencyMap = new AliMUONTriggerEfficiencyCells(filename, listname);
79   FillFromList();
80 }
81
82
83 //_____________________________________________________________________________
84 AliMUONTriggerChamberEfficiency::AliMUONTriggerChamberEfficiency(const AliMUONTriggerChamberEfficiency& other)
85 :
86 TObject(other),
87 fIsOwner(other.fIsOwner),
88 fEfficiencyMap(other.fEfficiencyMap),
89 fEfficiencyObjects(other.fEfficiencyObjects),
90 fDisplayList(other.fDisplayList)
91 {
92 /// Copy constructor
93 }
94
95 //_____________________________________________________________________________
96 AliMUONTriggerChamberEfficiency& AliMUONTriggerChamberEfficiency::operator=(const AliMUONTriggerChamberEfficiency& other)
97 {
98   /// Asignment operator
99   // check assignement to self
100   if (this == &other)
101     return *this;
102
103   fIsOwner = other.fIsOwner;
104   fEfficiencyMap = other.fEfficiencyMap;
105   fEfficiencyObjects = other.fEfficiencyObjects;
106   fDisplayList = other.fDisplayList;
107     
108   return *this;
109 }
110
111 //__________________________________________________________________________
112 AliMUONTriggerChamberEfficiency::~AliMUONTriggerChamberEfficiency()
113 {
114 ///  Destructor.
115   if ( fIsOwner )
116     delete fEfficiencyMap;
117   delete fEfficiencyObjects;
118   delete fDisplayList;
119 }
120
121
122 //__________________________________________________________________________
123 Float_t AliMUONTriggerChamberEfficiency::GetCellEfficiency(Int_t detElemId, Int_t localBoard, Int_t hType) const
124 {
125 ///  Get the efficiencies of the 2 cathodes at a given local board
126
127   Int_t chamber = FindChamberIndex(detElemId);
128   Int_t index = GetIndex(kHboardEff, hType, chamber);
129   TGraphAsymmErrors* effGraph = ((TGraphAsymmErrors*)fEfficiencyObjects->At(index));
130
131   // Some graphs are not available in the old implementation
132   if ( ! effGraph ) return -1.;
133
134   Double_t xpt, ypt;
135   effGraph->GetPoint(localBoard-1, xpt, ypt);
136   return ypt;
137 }
138
139
140 //__________________________________________________________________________
141 Float_t AliMUONTriggerChamberEfficiency::GetCellEfficiencyError(Int_t detElemId, Int_t localBoard, Int_t hType, Int_t errType) const
142 {
143   /// Get the efficiencie errors of the 2 cathodes at a given local board
144   /// errype 0 -> low error
145   /// errype 1 -> high error
146   
147   Int_t chamber = FindChamberIndex(detElemId);
148   Int_t index = GetIndex(kHboardEff, hType, chamber);
149   TGraphAsymmErrors* effGraph = ((TGraphAsymmErrors*)fEfficiencyObjects->At(index));
150   
151   // Some graphs are not available in the old implementation
152   if ( ! effGraph ) return -1.;
153
154   Float_t err = -1.;
155   if ( errType == 0 ) { // low error
156     err = effGraph->GetErrorYlow(localBoard-1);
157   }
158   else if ( errType == 1 ) { // up error
159     err = effGraph->GetErrorYhigh(localBoard-1);
160   }
161   
162   return err;
163 }
164
165
166 //__________________________________________________________________________
167 void 
168 AliMUONTriggerChamberEfficiency::IsTriggered(Int_t detElemId, Int_t localBoard, Bool_t &trigBend, Bool_t &trigNonBend) const
169 {
170 ///  Whether or not a given local board has a chance to trig, on each cathode.
171
172   // P(B) : probability to fire bending plane
173   Float_t effBend = GetCellEfficiency(detElemId, localBoard, AliMUONTriggerEfficiencyCells::kBendingEff);
174
175   // P(BN) : probability to fire bending and non-bending plane
176   Float_t effBoth = GetCellEfficiency(detElemId, localBoard, AliMUONTriggerEfficiencyCells::kBothPlanesEff);
177
178   trigBend =  ( gRandom->Rndm() > effBend ) ? kFALSE : kTRUE;
179
180   // P(N) : probability to fire non-bending plane
181   Float_t effNonBend = GetCellEfficiency(detElemId, localBoard, AliMUONTriggerEfficiencyCells::kNonBendingEff);
182
183   if ( effBoth > 0 ) {
184     effNonBend = ( trigBend ) ? 
185       effBoth / effBend :   // P(N|B) = P(BN) / P(B)
186       ( effNonBend - effBoth ) / ( 1. - effBend );  // P(N|!B) = ( P(N) - P(BN) ) / ( 1 - P(B) )
187   }
188
189   trigNonBend =  ( gRandom->Rndm() > effNonBend ) ? kFALSE : kTRUE;
190
191   AliDebug(2,Form("Ch %i  board %i  resp (%i, %i)  prob (%.2f, %.2f)  effNB %.2f  effBoth %.2f\n", detElemId/100, localBoard, trigBend, trigNonBend, effBend, effNonBend, GetCellEfficiency(detElemId, localBoard, AliMUONTriggerEfficiencyCells::kNonBendingEff), effBoth));
192
193
194 }
195
196
197 //__________________________________________________________________________
198 Int_t AliMUONTriggerChamberEfficiency::FindChamberIndex(Int_t detElemId) const
199 {
200 ///  From detElemId to chamber number
201
202   // Int_t iChamber = AliMpDEManager::GetChamberId(detElemId);
203   Int_t iChamber = detElemId/100 - 1;
204   return iChamber-AliMpConstants::NofTrackingChambers();
205 }
206
207
208 //__________________________________________________________________________
209 void
210 AliMUONTriggerChamberEfficiency::FillFromList(Bool_t useMeanValues)
211 {
212 ///  Fills internal histos from list.
213
214   if ( fEfficiencyObjects )
215     delete fEfficiencyObjects;
216
217   const Int_t kNeffHistos = 
218     2 * ( AliMUONTriggerEfficiencyCells::kNcounts - 1 ) * AliMpConstants::NofTriggerChambers();
219
220   fEfficiencyObjects = new TObjArray(kNeffHistos);
221   fEfficiencyObjects->SetOwner();
222
223   TH1F *histoNum = 0x0, *histoDen=0x0;
224   TString histoName = "";
225   Int_t deType[2] = {AliMUONTriggerEfficiencyCells::kHboardCount,
226                      AliMUONTriggerEfficiencyCells::kHslatCount};
227   Int_t deTypeEff[2] = {kHboardEff, kHslatEff};
228   Int_t index = -1;
229
230   Bool_t rebuildEfficiency = kTRUE;
231
232   for ( Int_t ide=0; ide<2; ide++){
233     Int_t currDe = deType[ide];
234
235     if ( useMeanValues && currDe == AliMUONTriggerEfficiencyCells::kHboardCount ) 
236       continue;
237
238     for(Int_t ich=0; ich<AliMpConstants::NofTriggerChambers(); ich++){
239       histoName = fEfficiencyMap->GetHistoName(currDe, AliMUONTriggerEfficiencyCells::kAllTracks, ich);
240       if ( fEfficiencyMap->GetHistoList() ) {
241         histoDen = (TH1F*)fEfficiencyMap->GetHistoList()->FindObject(histoName.Data());
242         if ( !histoDen ) {
243           AliWarning(Form("Histogram %s not found. Efficiency won't be re-build", histoName.Data()));
244           rebuildEfficiency = kFALSE;
245         }
246       }
247       else {
248         AliWarning("Histogram list not present: efficiency won't be re-build");
249         rebuildEfficiency = kFALSE;
250       }
251
252       Int_t nTypes = ( rebuildEfficiency ) ? AliMUONTriggerEfficiencyCells::kNcounts-1 : 2; 
253       for(Int_t hType=0; hType<nTypes; hType++){
254         histoName = fEfficiencyMap->GetHistoName(currDe, hType, ich);
255
256         histoNum = ( rebuildEfficiency ) ? 
257           (TH1F*)fEfficiencyMap->GetHistoList()->FindObject(histoName.Data()) :
258           fEfficiencyMap->GetOldEffHisto(currDe, ich, hType);
259       
260         if ( !histoNum ) {
261           AliWarning(Form("Histogram %s not found. Skip to next", histoName.Data()));
262           continue;
263         }
264
265         index = GetIndex(deTypeEff[ide], hType, ich);
266         TGraphAsymmErrors* effGraph = GetEfficiencyGraph(histoNum,histoDen);
267         fEfficiencyObjects->AddAt(effGraph, index);
268
269         TString debugString = Form("Adding object %s",effGraph->GetName());
270         if ( histoDen ) debugString += Form(" (%s/%s)",histoNum->GetName(),histoDen->GetName());
271         debugString += Form(" index %i",index);
272         AliDebug(5,debugString.Data());
273
274         if ( useMeanValues && rebuildEfficiency ){
275           Int_t currChamber = ich + AliMpConstants::NofTrackingChambers();
276           histoName = fEfficiencyMap->GetHistoName(AliMUONTriggerEfficiencyCells::kHboardCount, hType, ich);
277           TH1F* auxHistoNum = (TH1F*)fEfficiencyMap->GetHistoList()->FindObject(histoName.Data())->Clone("tempHistoNum");
278           TH1F* auxHistoDen = (TH1F*)fEfficiencyMap->GetHistoList()->FindObject(histoName.Data())->Clone("tempHistoDen");
279           for ( Int_t iBinBoard = 1; iBinBoard<=AliMpConstants::NofLocalBoards(); iBinBoard++){
280             Int_t detElemId = AliMpDDLStore::Instance()->GetDEfromLocalBoard(iBinBoard, currChamber);
281             Int_t iBin = histoNum->FindBin(detElemId%100);
282
283             auxHistoNum->SetBinContent(iBinBoard, histoNum->GetBinContent(iBin));
284             auxHistoDen->SetBinContent(iBinBoard, histoDen->GetBinContent(iBin));
285           }
286           index = GetIndex(kHboardEff, hType, ich);
287           effGraph = GetEfficiencyGraph(auxHistoNum,auxHistoDen);
288           fEfficiencyObjects->AddAt(effGraph, index);
289           AliDebug(5,Form("Adding object %s (%s/%s) at index %i",effGraph->GetName(),histoNum->GetName(),histoDen->GetName(),index));
290           delete auxHistoNum;
291           delete auxHistoDen;
292         } // if (useMeanValues)
293       } // loop on count type
294     } // loop on chamber
295   } // loop on detection element histogram
296 }
297
298
299 //_____________________________________________________________________________
300 void AliMUONTriggerChamberEfficiency::DisplayEfficiency(Bool_t perSlat, Bool_t show2Dhisto)
301 {
302   //
303   /// Display calculated efficiency.
304   //
305
306   if ( !AliCDBManager::Instance()->GetDefaultStorage() ){
307     AliWarning("Please set default CDB storage (needed for mapping).");
308     return;
309   }
310   if ( AliCDBManager::Instance()->GetRun() < 0 ){
311     AliWarning("Please set CDB run number (needed for mapping).");
312     return;
313   }
314
315   TString baseCanName = "MTRtrigChEffCan";
316   TString histoName;
317
318   // Remove previously created canvases
319   TCanvas* can = 0x0;
320   TIter next(gROOT->GetListOfCanvases());
321   while ((can = (TCanvas *)next())) {
322     histoName = can->GetName();
323     if ( histoName.Contains(baseCanName.Data()))
324       delete can;
325   }
326
327   delete fDisplayList;
328   fDisplayList = new TList();
329   fDisplayList->SetOwner();
330
331   TH2F* displayHisto = 0x0;
332
333   AliMUONTriggerDisplay triggerDisplay;
334
335   Int_t deType = ( perSlat ) ? kHslatEff : kHboardEff;
336   AliMUONTriggerDisplay::EDisplayType displayType = ( perSlat ) ? 
337     AliMUONTriggerDisplay::kDisplaySlats : AliMUONTriggerDisplay::kDisplayBoards;
338   Int_t index = -1;
339
340   TGraph* graph = 0x0;
341
342   // Book histos
343   for(Int_t ich=0; ich<AliMpConstants::NofTriggerChambers(); ich++){
344     Int_t currCh = 11 + ich;
345     for(Int_t hType=0; hType<AliMUONTriggerEfficiencyCells::kNcounts - 1; hType++){
346       index = GetIndex(deType, hType, ich);
347       graph = (TGraph*)fEfficiencyObjects->At(index);
348       if ( ! graph ) continue;
349       histoName = graph->GetName();
350       histoName += baseCanName;
351       Int_t shift = 10*(index%((AliMUONTriggerEfficiencyCells::kNcounts - 1)*
352                                AliMpConstants::NofTriggerChambers()));
353       can = new TCanvas(histoName.Data(), histoName.Data(), 100+shift, shift, 700, 700);
354       can->SetRightMargin(0.14);
355       can->SetLeftMargin(0.12);
356       histoName.ReplaceAll(baseCanName.Data(), "Display");
357       if ( show2Dhisto ) {
358         displayHisto = 
359           (TH2F*)triggerDisplay.GetDisplayHistogram(graph, histoName,
360                                                     displayType,
361                                                     hType,currCh,histoName,
362                                                     AliMUONTriggerDisplay::kShowZeroes);
363         displayHisto->SetDirectory(0);
364       }
365
366       if ( show2Dhisto ){
367         displayHisto->GetZaxis()->SetRangeUser(0.,1.);
368         displayHisto->GetYaxis()->SetTitleOffset(1.4);
369         displayHisto->SetStats(kFALSE);
370         displayHisto->DrawCopy("COLZ");
371         delete displayHisto;
372
373         if ( deType == kHboardEff ){
374           histoName = Form("labels%iChamber%i", hType, currCh);
375           displayHisto = 
376             (TH2F*)triggerDisplay.GetBoardNumberHisto(histoName,currCh);
377           displayHisto->SetDirectory(0);
378           displayHisto->DrawCopy("textsame");
379           delete displayHisto;
380         }
381       }
382       else {
383         TGraphAsymmErrors* drawGraph = (TGraphAsymmErrors*)graph->Clone(histoName.Data());
384         drawGraph->SetMarkerStyle(20);
385         drawGraph->SetMarkerSize(0.7);
386         drawGraph->SetMarkerColor(kRed);
387         fDisplayList->Add(drawGraph);
388         drawGraph->Draw("ap");
389       } // loop on chamber
390     } // loop on count type
391   } // loop on chamber
392 }
393
394
395 //__________________________________________________________________________
396 Bool_t AliMUONTriggerChamberEfficiency::LowStatisticsSettings(Bool_t useMeanValues)
397 {
398   //
399   /// In case of low statistics, fill the local board efficiency with
400   /// the average value of the RPC
401   //
402
403   if ( useMeanValues )
404     AliInfo("Boards filled with the average efficiency of the RPC");
405
406   FillFromList(useMeanValues);
407
408   return kTRUE;
409 }
410
411
412 //__________________________________________________________________________
413 Int_t
414 AliMUONTriggerChamberEfficiency::GetIndex(Int_t histoType, Int_t countType, 
415                                           Int_t chamber) const
416 {
417   //
418   /// Return the index of the object in the array
419   //
420
421   const Int_t kNtypes = AliMUONTriggerEfficiencyCells::kNcounts - 1;
422   const Int_t kNchambers = AliMpConstants::NofTriggerChambers();
423   return 
424     histoType * kNtypes * kNchambers + 
425     chamber * kNtypes +
426     countType;
427   //countType * kNchambers +
428   //chamber;
429 }
430
431
432 //_____________________________________________________________________________
433 TObject* AliMUONTriggerChamberEfficiency::GetEffObject(Int_t histoType, Int_t countType, 
434                                                                  Int_t chamber)
435 {
436   //
437   /// Get efficiency object
438   //
439   Int_t index = GetIndex(histoType, countType, chamber);
440   return fEfficiencyObjects->At(index);
441 }
442
443
444 //_____________________________________________________________________________
445 TGraphAsymmErrors* AliMUONTriggerChamberEfficiency::GetEfficiencyGraph(TH1* histoNum, TH1* histoDen)
446 {
447   //
448   /// Create the graph of efficiency from the numerator and denominator
449   /// histogram in such a way to have a point set also for
450   /// detection elements with efficiency = 0 or non calculated
451   //
452
453   TGraphAsymmErrors* auxGraph = 0x0;
454   if ( histoDen ) auxGraph = new TGraphAsymmErrors(histoNum,histoDen,"cp");
455   else auxGraph = new TGraphAsymmErrors(histoNum);
456
457   Int_t npoints = histoNum->GetNbinsX();
458   TGraphAsymmErrors* effGraph = new TGraphAsymmErrors(npoints);
459   TString histoName = histoNum->GetName();
460   histoName.ReplaceAll("Count","Eff");
461   effGraph->SetName(histoName.Data());
462   effGraph->SetTitle(histoName.Data());
463   Double_t oldX, oldY;
464   for ( Int_t ibin=0; ibin<npoints; ibin++ ) {
465     Int_t foundPoint = -1;
466     for (Int_t ipt=0; ipt<auxGraph->GetN(); ipt++) {
467       auxGraph->GetPoint(ipt, oldX, oldY);
468       if ( oldX > histoNum->GetBinLowEdge(ibin+1) &&
469            oldX < histoNum->GetBinLowEdge(ibin+2) ) {
470         foundPoint = ipt;
471         break;
472       }
473     }
474     Double_t currX = ( foundPoint < 0 ) ? histoNum->GetBinCenter(ibin+1) : oldX; 
475     Double_t currY = ( foundPoint < 0 ) ? 0. : oldY;
476     Double_t eyl   = ( foundPoint < 0 ) ? 0. : auxGraph->GetErrorYlow(foundPoint);
477     Double_t eyh   = ( foundPoint < 0 ) ? 0. : auxGraph->GetErrorYhigh(foundPoint);
478     effGraph->SetPoint(ibin, currX, currY);
479     effGraph->SetPointError(ibin, 0., 0., eyl, eyh);
480   }
481
482   delete auxGraph;
483
484   return effGraph;
485 }