827c65c0616ffa697f82cda7c5ac931343dcb4f1
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackerQAChecker.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 "AliMUONTrackerQAChecker.h"
19
20 /// \class AliMUONTrackerQAChecker
21 ///
22 /// Implementation of AliQACheckerBase for MCH and MTR
23 ///
24 /// For the moment we only implement the checking of raw data QA for the tracker
25 /// by looking at the occupancy at the bus patch level.
26 ///
27 /// \author Laurent Aphecetche, Subatech
28
29 #include "AliCDBManager.h"
30 #include "AliCodeTimer.h"
31 #include "AliLog.h"
32 #include "AliMUONQAIndices.h"
33 #include "AliMUONRecoParam.h"
34 #include "AliMpBusPatch.h"
35 #include "AliMpDDLStore.h"
36 #include "AliQAv1.h"
37 #include "AliQAv1.h"
38 #include "Riostream.h"
39 #include "TAxis.h"
40 #include "TDirectory.h"
41 #include "TH1.h"
42 #include "TLine.h"
43 #include "TMath.h"
44 #include "TPaveText.h"
45
46 /// \cond CLASSIMP
47 ClassImp(AliMUONTrackerQAChecker)
48 /// \endcond
49
50 namespace {
51   
52   //___________________________________________________________________________  
53   int trim(Int_t n, 
54            Double_t* x,
55            Double_t alpha,
56            Double_t& tmean,
57            Double_t& tvar,
58            Double_t& min,
59            Double_t& max)
60   {
61     //
62     // Calculates the trimmed (tmean) mean
63     // of a sample (x) and estimates the variance (tvar)
64     // of that mean.
65     //
66     
67     // First check input parameters
68     
69     // number of observations
70     if ( n < 2 )
71     {
72       return -1;
73     }
74     
75     if ( alpha < 0 || alpha >= 0.5 )
76       // proportion of observations
77       // to be trimmed at each end of the sorted sample
78     {
79       return -2;
80     }
81     
82     // Input parameters are good. Let's move on.
83     
84     // Insure we use a sample sorted into ascending order.
85     
86     Int_t* indices = new Int_t[n];
87     
88     TMath::Sort(n,x,indices,kFALSE);
89     
90     Double_t* sx = new Double_t[n];
91     
92     for ( Int_t i = 0; i < n; ++i )
93     {
94       sx[i] = x[indices[i]];
95     }
96     delete[] indices;
97     
98     
99     // Number of observations trimmed at each end.
100     
101     Int_t k = TMath::FloorNint(alpha * n);
102     
103     double sum = 0.0;
104     
105     for ( Int_t i = k; i < n - k ; ++i )
106     {
107       sum += sx[i];
108     }
109     
110     tmean = sum / ( n - 2 * k );
111     
112     double t2 = 0.0;
113     
114     for ( Int_t i = k; i < n - k; ++i )
115     {
116       t2 += (sx[i] - tmean) * (sx[i] - tmean);
117     }
118     
119     tvar = (
120             t2 +
121             k * (sx[k] - tmean) * (sx[k] - tmean) +
122             k * (sx[n - k - 1] - tmean) * (sx[n - k - 1] - tmean)
123             ) / (n * n);
124     
125     // get the min and max for the non-rejected values
126     min = DBL_MAX;
127     max = 0.0;
128     
129     for ( Int_t i = k; i < n-k; ++i ) 
130     {
131       min = TMath::Min(min,sx[i]);
132       max = TMath::Max(max,sx[i]);
133     }
134     
135     delete[] sx;
136     
137     return 0;
138   }
139 }
140
141 //__________________________________________________________________
142 AliMUONTrackerQAChecker::AliMUONTrackerQAChecker() : AliMUONVQAChecker()
143 {
144         /// ctor
145 }          
146
147 //__________________________________________________________________
148 AliMUONTrackerQAChecker::~AliMUONTrackerQAChecker() 
149 {
150         /// dtor
151 }
152
153 //__________________________________________________________________
154 AliMUONTrackerQAChecker::AliMUONTrackerQAChecker(const AliMUONTrackerQAChecker& qac) : 
155     AliMUONVQAChecker(qac) 
156 {
157         /// copy ctor 
158 }   
159
160 //______________________________________________________________________________
161 AliMUONVQAChecker::ECheckCode*
162 AliMUONTrackerQAChecker::CheckRecPoints(TObjArray ** list, const AliMUONRecoParam* /*recoParam*/)
163 {
164   /// Check rec points
165   /// Very binary check for the moment. 
166   
167   AliCodeTimerAuto("",0);
168   
169   AliMUONVQAChecker::ECheckCode * rv = new AliMUONVQAChecker::ECheckCode[AliRecoParam::kNSpecies] ; 
170   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
171     rv[specie] =  AliMUONVQAChecker::kInfo; 
172   
173   
174   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
175   {
176     TH1* h = AliQAv1::GetData(list,AliMUONQAIndices::kTrackerNumberOfClustersPerDE,AliRecoParam::ConvertIndex(specie));
177
178     if ( !h ) rv[specie] = AliMUONVQAChecker::kWarning; // only a warning if histo not found, in order not to kill anything because QA failed...
179   
180     else if ( h->GetMean() == 0.0 ) rv[specie] =  MarkHisto(*h,AliMUONVQAChecker::kFatal);
181   }
182   return rv;
183 }
184
185 //______________________________________________________________________________
186 AliMUONVQAChecker::ECheckCode 
187 AliMUONTrackerQAChecker::MarkHisto(TH1& histo, AliMUONVQAChecker::ECheckCode value) const
188 {
189   /// Mark histo as originator of some QA error/warning
190   
191   if ( value != AliMUONVQAChecker::kInfo )
192   {
193     histo.SetBit(AliQAv1::GetQABit());
194   }
195   
196   return value;
197 }
198
199 //______________________________________________________________________________
200 AliMUONVQAChecker::ECheckCode*
201 AliMUONTrackerQAChecker::CheckESD(TObjArray ** list, const AliMUONRecoParam* /*recoParam*/)
202 {
203   /// Check ESD
204   
205   AliCodeTimerAuto("",0);
206   
207   AliMUONVQAChecker::ECheckCode * rv = new AliMUONVQAChecker::ECheckCode[AliRecoParam::kNSpecies] ; 
208   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
209     rv[specie] = AliMUONVQAChecker::kInfo; 
210   
211   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
212     
213     TH1* h = AliQAv1::GetData(list,AliMUONQAIndices::kESDnTracks,AliRecoParam::ConvertIndex(specie));
214   
215     if (!h) rv[specie] = AliMUONVQAChecker::kWarning;
216   
217     else if ( h->GetMean() == 0.0 ) rv[specie] =  MarkHisto(*h,AliMUONVQAChecker::kFatal); // no track -> fatal
218   
219     h = AliQAv1::GetData(list,AliMUONQAIndices::kESDMatchTrig,AliRecoParam::ConvertIndex(specie));
220   
221     if (!h) rv[specie] = AliMUONVQAChecker::kWarning;
222   
223     else if (h->GetMean() == 0.0 ) rv[specie] = MarkHisto(*h,AliMUONVQAChecker::kError); // no trigger matching -> error
224   }
225   return rv;
226 }
227
228 //______________________________________________________________________________
229 AliMUONVQAChecker::ECheckCode*
230 AliMUONTrackerQAChecker::CheckRaws(TObjArray ** list, const AliMUONRecoParam* recoParam)
231 {
232   /// Check raws
233
234   AliCodeTimerAuto("",0);
235   
236   if (!recoParam) return 0x0;
237   
238   AliMUONVQAChecker::ECheckCode * rv = new AliMUONVQAChecker::ECheckCode[AliRecoParam::kNSpecies] ; 
239
240   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
241   {
242     rv[specie] = AliMUONVQAChecker::kInfo; 
243   }
244   
245   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
246   {
247     TH1* hbp = AliQAv1::GetData(list,AliMUONQAIndices::kTrackerBusPatchOccupancy,AliRecoParam::ConvertIndex(specie));    
248
249     TH1* hnpads = AliQAv1::GetData(list,AliMUONQAIndices::kTrackerBusPatchNofPads,AliRecoParam::ConvertIndex(specie));
250
251     TH1* hbpconfig = AliQAv1::GetData(list,AliMUONQAIndices::kTrackerBusPatchConfig,AliRecoParam::ConvertIndex(specie));
252
253     TH1* hnevents = AliQAv1::GetData(list,AliMUONQAIndices::kTrackerNofRawEventSeen,AliRecoParam::ConvertIndex(specie));
254
255     if ( !hbp || !hnpads || !hnevents ) 
256     {
257       continue;
258     }
259
260     Int_t nevents = TMath::Nint(hnevents->GetBinContent(1));
261     
262     rv[specie] = BeautifyTrackerBusPatchOccupancy(*hbp,hbpconfig,*hnpads,nevents,*recoParam);    
263   }
264   
265   return rv;
266 }
267
268 //____________________________________________________________________________ 
269 AliMUONVQAChecker::ECheckCode
270 AliMUONTrackerQAChecker::BeautifyTrackerBusPatchOccupancy(TH1& hbp, 
271                                                           const TH1* hbuspatchconfig, 
272                                                           const TH1& hnpads, 
273                                                           Int_t nevents,
274                                                           const AliMUONRecoParam& recoParam)
275 {
276   /// Put labels, limits and so on on the TrackerBusPatchOccupancy histogram
277   /// hbuspatchconfig and hbp must have the same bin definitions
278   
279   if ( hbuspatchconfig )
280   {
281     if ( hbp.GetNbinsX() != hbuspatchconfig->GetNbinsX() ||
282         hbp.GetXaxis()->GetXmin() != hbuspatchconfig->GetXaxis()->GetXmin() ||
283         hbp.GetXaxis()->GetXmax() != hbuspatchconfig->GetXaxis()->GetXmax() )
284     {
285       AliError("hbp and hbuspatchconfig histograms are not compatible !");
286       return AliMUONVQAChecker::kFatal;
287     }
288   }
289   
290   hbp.SetXTitle("Absolute Bus Patch Id");
291   hbp.SetYTitle("Occupancy (percent)");
292   hbp.SetStats(kFALSE);
293   
294   Double_t xmin = hbp.GetXaxis()->GetXmin();
295   Double_t xmax = hbp.GetXaxis()->GetXmax();
296   
297   Double_t occMax(0.1); // 0.1% y-limit for the plot
298   Double_t maxToleratedOccupancy(recoParam.BuspatchOccupancyHighLimit()*100.0); 
299   Double_t minToleratedOccupancy(recoParam.BuspatchOccupancyLowLimit()*100.0);   
300   TLine* line1 = new TLine(xmin,maxToleratedOccupancy,xmax,maxToleratedOccupancy);
301   line1->SetLineColor(1);
302   line1->SetLineWidth(1);
303
304   TLine* line2 = new TLine(xmin,minToleratedOccupancy,xmax,minToleratedOccupancy);
305   line2->SetLineColor(1);
306   line2->SetLineWidth(1);
307   
308   hbp.GetListOfFunctions()->Add(line1);
309   hbp.GetListOfFunctions()->Add(line2);
310     
311   TIter next(AliMpDDLStore::Instance()->CreateBusPatchIterator());
312   AliMpBusPatch* bp(0x0);
313
314   Int_t nMissingPads(0);
315   Int_t nPads(0);
316   Int_t nBusPatches(0);
317   Int_t nMissingBusPatches(0);
318   
319   while ( ( bp = static_cast<AliMpBusPatch*>(next())) )
320   {
321     Int_t bin = hbp.FindBin(bp->GetId());
322     Int_t n = TMath::Nint(hnpads.GetBinContent(bin));
323     
324     ++nBusPatches;
325     
326     nPads += n;
327     
328     if ( hbp.GetBinContent(bin) <= 0 ) 
329     {
330       nMissingPads += n;
331       ++nMissingBusPatches;
332     }
333   }
334   
335   next.Reset();
336   
337   Int_t ok(-1);
338   Int_t n(0);
339   Int_t nBusPatchesAboveLimit(0);
340   Int_t nBusPatchesBelowLimit(0);
341   Double_t alpha(0.1); // trim 10% of data
342   Double_t tmean(0.0),tvar(0.0);
343   Double_t ymin(0.0),ymax(0.0);
344   AliMUONVQAChecker::ECheckCode rv(AliMUONVQAChecker::kFatal); // default value = serious problem
345   
346   if ( nBusPatches ) 
347   {
348     Double_t* x = new Double_t[nBusPatches];
349     
350     while ( ( bp = static_cast<AliMpBusPatch*>(next())) )
351     {
352       Int_t bin = hbp.FindBin(bp->GetId());
353       if ( hbp.GetBinContent(bin) > 0 )
354       {
355         x[n] = hbp.GetBinContent(bin);
356         ++n;
357       }
358       if ( hbp.GetBinContent(bin) > maxToleratedOccupancy )
359       {
360         ++nBusPatchesAboveLimit;
361       }
362       if ( hbp.GetBinContent(bin) < minToleratedOccupancy )
363       {
364         // check whether this buspatch has a reason to be absent (only valid
365         // if we got the config, otherwise we cannot do the test)
366         if ( hbuspatchconfig && hbuspatchconfig->GetBinContent(bin) > 0 )
367         {
368           // should be there, so it's an error
369           ++nBusPatchesBelowLimit;
370         }
371       }
372     }
373     
374     // computed the truncated mean of the occupancy values, in order to get a 
375     // reasonable y-range for the histogram (without giant peaks to the roof 
376     // for misbehaving buspatches).
377     ok = trim(n,x,alpha,tmean,tvar,ymin,ymax);
378     
379     delete[] x;
380   }
381   
382   if ( ok < 0 ) 
383   {
384     ymax = occMax;
385   }
386   else
387   {
388     ymax = TMath::Max(ymax,occMax);
389   }
390   
391   hbp.SetMaximum(ymax*1.4);
392   
393   TPaveText* text = new TPaveText(0.50,0.60,0.99,0.99,"NDC");
394   
395   text->AddText(Form("MCH RUN %d - %d events",AliCDBManager::Instance()->GetRun(),nevents));
396   
397   if ( ok < 0 ) 
398   {
399     text->AddText("Could not compute truncated mean. Not enough events ?");
400     text->AddText(Form("nBusPatches=%d n=%d",nBusPatches,n));
401   }
402   else if (!nPads || !nBusPatches)
403   {
404     text->AddText(Form("Could not get the total number of pads (%d) or total number of buspatches (%d). ERROR !!!",
405                        nPads,nBusPatches));
406   }
407   else
408   {
409     Float_t missingPadFraction = nMissingPads*100.0/nPads;
410     Float_t missingBusPatchFraction = nMissingBusPatches*100.0/nBusPatches;
411     Float_t aboveLimitFraction = nBusPatchesAboveLimit*100.0/nBusPatches;
412     Float_t belowLimitFraction = nBusPatchesBelowLimit*100.0/nBusPatches;
413     
414     text->AddText(Form("%5.2f %% of missing buspatches (%d out of %d)",missingBusPatchFraction,nMissingBusPatches,nBusPatches));
415     text->AddText(Form("%5.2f %% of missing pads (%d out of %d)",missingPadFraction,nMissingPads,nPads));
416     text->AddText(Form("%5.2f %% bus patches above the %5.2f %% limit",aboveLimitFraction,maxToleratedOccupancy));
417     text->AddText(Form("%5.2f %% bus patches below the %e %% limit",belowLimitFraction,minToleratedOccupancy));
418     text->AddText(Form("Truncated mean at %2d %% is %7.2f %%",(Int_t)(alpha*100),tmean));
419     
420     if ( missingPadFraction >= 100.0 ) 
421     {
422       rv = AliMUONVQAChecker::kFatal;       
423     }
424     
425     else if ( missingPadFraction > recoParam.MissingPadFractionLimit()*100.0 || 
426              aboveLimitFraction > recoParam.FractionOfBuspatchOutsideOccupancyLimit()*100.0 || 
427              belowLimitFraction > recoParam.FractionOfBuspatchOutsideOccupancyLimit()*100.0 ) 
428     {
429       rv = AliMUONVQAChecker::kError;
430     }
431     else
432     {
433       rv = AliMUONVQAChecker::kInfo;
434     }
435   }
436   
437   hbp.GetListOfFunctions()->Add(text);
438   
439   if ( rv == AliMUONVQAChecker::kInfo ) 
440   {
441     text->SetFillColor(3); // green = INFO
442   }
443   else if ( rv == AliMUONVQAChecker::kWarning ) 
444   {
445     text->SetFillColor(5); // yellow = WARNING
446   }
447   else if ( rv ==  AliMUONVQAChecker::kFatal) 
448   {
449     text->SetFillColor(2); // red = FATAL
450   }
451   else 
452   {
453     text->SetFillColor(6); // pink = ERROR
454   }
455   
456   return rv;
457 }
458