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