Protection against division by zero.
[u/mrichter/AliRoot.git] / TPC / AliTPCQADataMakerRec.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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
17 /* $Id: $ */
18
19 /*
20   Based on AliPHOSQADataMaker
21   Produces the data needed to calculate the quality assurance. 
22   All data must be mergeable objects.
23   P. Christiansen, Lund, January 2008
24 */
25
26 /*
27   Implementation:
28
29   We have chosen to have the histograms as non-persistent meber to
30   allow better debugging. In the copy constructor we then have to
31   assign the pointers to the existing histograms in the copied
32   list. This have been implemented but not tested.
33
34   For the QA of the RAW data we use the class, AliTPCdataQA, from the
35   existing TPC Calibration framework (which is more advanced than the
36   standard QA framework) and extract the histograms at the end. This
37   has been tested with zero-suppressed data. The Analyse method of the
38   AliTPCdataQA class is called in the method, EndOfDetectorCycle, and
39   there also: 1d histogram(s) are projected and added to the QA list.
40 */
41
42 /*
43   TODO:
44   Sumw2 for RAW histogram(s)?
45   RecPoints and ESD could have many more histograms
46 */
47
48 #include "AliTPCQADataMakerRec.h"
49
50 // --- ROOT system ---
51 #include <TClonesArray.h>
52 #include <TString.h>
53 #include <TSystem.h>
54
55 // --- Standard library ---
56
57 // --- AliRoot header files ---
58 #include "AliQAChecker.h"
59 #include "AliESDEvent.h"
60 #include "AliESDtrack.h"
61 #include "AliLog.h"
62 #include "AliTPCCalPad.h"
63 #include "AliTPCCalROC.h"
64 #include "AliTPCClustersRow.h"
65 #include "AliTPCclusterMI.h"
66
67 ClassImp(AliTPCQADataMakerRec)
68
69 //____________________________________________________________________________ 
70 AliTPCQADataMakerRec::AliTPCQADataMakerRec() : 
71   AliQADataMakerRec(AliQA::GetDetName(AliQA::kTPC), 
72                     "TPC Rec Quality Assurance Data Maker"),
73   fTPCdataQA(0),
74   fHistESDclusters(0),fHistESDratio(0), fHistESDpt(0),
75   fHistRawsOccupancy(0), fHistRawsOccupancyVsSector(0),
76   fHistRawsNClustersPerEventVsSector(0), fHistRawsQVsSector(0),
77   fHistRawsQmaxVsSector(0),
78   fHistRecPointsQmaxShort(0), fHistRecPointsQmaxMedium(0), 
79   fHistRecPointsQmaxLong(0), fHistRecPointsQShort(0), 
80   fHistRecPointsQMedium(0), fHistRecPointsQLong(0),
81   fHistRecPointsRow(0)
82 {
83   // ctor
84   for(Int_t i = 0; i < 6; i++)
85     fMapping[i] = 0;
86 }
87
88 //____________________________________________________________________________ 
89 AliTPCQADataMakerRec::AliTPCQADataMakerRec(const AliTPCQADataMakerRec& qadm) :
90   AliQADataMakerRec(),
91   fTPCdataQA(0),
92   fHistESDclusters(0),fHistESDratio(0), fHistESDpt(0),
93   fHistRawsOccupancy(0), fHistRawsOccupancyVsSector(0),
94   fHistRawsNClustersPerEventVsSector(0), fHistRawsQVsSector(0),
95   fHistRawsQmaxVsSector(0),
96   fHistRecPointsQmaxShort(0), fHistRecPointsQmaxMedium(0), 
97   fHistRecPointsQmaxLong(0), fHistRecPointsQShort(0), 
98   fHistRecPointsQMedium(0), fHistRecPointsQLong(0),
99   fHistRecPointsRow(0)
100 {
101   //copy ctor 
102   // Does not copy the calibration object, instead InitRaws have to be
103   // called again
104   SetName((const char*)qadm.GetName()) ; 
105   SetTitle((const char*)qadm.GetTitle()); 
106
107   for(Int_t i = 0; i < 6; i++)
108     fMapping[i] = 0;
109
110   //
111   // Associate class histogram objects to the copies in the list
112   // Could also be done with the indexes
113   //
114   fHistESDclusters   = (TH1F*)fESDsQAList->FindObject("hESDclusters");
115   fHistESDratio      = (TH1F*)fESDsQAList->FindObject("hESDratio");
116   fHistESDpt         = (TH1F*)fESDsQAList->FindObject("hESDpt");
117
118   fHistRawsOccupancy = (TH1F*)fRawsQAList->FindObject("hRawsOccupancy");
119   fHistRawsOccupancyVsSector = 
120     (TH1F*)fRawsQAList->FindObject("hRawsOccupancyVsSector");
121   fHistRawsNClustersPerEventVsSector = 
122     (TH1F*)fRawsQAList->FindObject("hRawsNClustersPerEventVsSector");
123   fHistRawsQVsSector = (TH1F*)fRawsQAList->FindObject("hRawsQVsSector");
124   fHistRawsQmaxVsSector = (TH1F*)fRawsQAList->FindObject("hRawsQmaxVsSector");
125   
126   fHistRecPointsQmaxShort  = 
127     (TH1F*)fRecPointsQAList->FindObject("hRecPointsQmaxShort");
128   fHistRecPointsQmaxMedium = 
129     (TH1F*)fRecPointsQAList->FindObject("hRecPointsQmaxMedium");
130   fHistRecPointsQmaxLong   = 
131     (TH1F*)fRecPointsQAList->FindObject("hRecPointsQmaxLong"); 
132   fHistRecPointsQShort     = 
133     (TH1F*)fRecPointsQAList->FindObject("hRecPointsQShort");
134   fHistRecPointsQMedium    = 
135     (TH1F*)fRecPointsQAList->FindObject("hRecPointsQMedium");
136   fHistRecPointsQLong      = 
137     (TH1F*)fRecPointsQAList->FindObject("hRecPointsQLong");
138   fHistRecPointsRow      = 
139     (TH1F*)fRecPointsQAList->FindObject("hRecPointsRow");
140 }
141
142 //__________________________________________________________________
143 AliTPCQADataMakerRec& AliTPCQADataMakerRec::operator = (const AliTPCQADataMakerRec& qadm )
144 {
145   // Equal operator.
146   this->~AliTPCQADataMakerRec();
147   new(this) AliTPCQADataMakerRec(qadm);
148   return *this;
149 }
150
151 //__________________________________________________________________
152 AliTPCQADataMakerRec::~AliTPCQADataMakerRec()
153 {
154   // Destructor
155   delete fTPCdataQA; 
156
157   for(Int_t i = 0; i < 6; i++) 
158     delete fMapping[i];
159 }
160  
161 //____________________________________________________________________________ 
162 void AliTPCQADataMakerRec::EndOfDetectorCycle(AliQA::TASKINDEX_t task, TObjArray * list)
163 {
164   //Detector specific actions at end of cycle
165
166   if(fTPCdataQA) { // do the final step of the QA for Raw data
167
168     fTPCdataQA->Analyse(); // 31/1-08 Analyse is now protected against
169                            //         RAW data files with no TPC data
170     
171     //Add2RawsList(fTPCdataQA, 0);
172     // get the histograms and add them to the output
173     // 31/8-08 Histogram is only added if the Calibration class 
174     //         receives TPC data 
175     const Int_t eventCounter = fTPCdataQA->GetEventCounter();
176     if(eventCounter>0) { // some TPC data has been processed
177
178       // Reset histograms and refill them 
179       fHistRawsOccupancy->Reset();
180       fHistRawsOccupancyVsSector->Reset();
181       fHistRawsNClustersPerEventVsSector->Reset();
182       fHistRawsQVsSector->Reset();
183       fHistRawsQmaxVsSector->Reset();
184       
185       TH1F* hNorm72 = new TH1F("hNorm72", "histogram to normalize 72 sectors",
186                                72, 0, 72);
187       hNorm72->Sumw2();
188       TH1F* hNorm108 = new TH1F("hNorm108", "histogram to normalize 108 sectors (medium and long pads are split up)",
189                                108, 0, 108);
190       hNorm108->Sumw2();
191
192       for (Int_t iSec = 0; iSec < 72; iSec++) {
193         
194         AliTPCCalROC* occupancyROC = 
195           fTPCdataQA->GetNoThreshold()->GetCalROC(iSec); 
196         AliTPCCalROC* nclusterROC = 
197           fTPCdataQA->GetNLocalMaxima()->GetCalROC(iSec); 
198         AliTPCCalROC* qROC = 
199           fTPCdataQA->GetMeanCharge()->GetCalROC(iSec); 
200         AliTPCCalROC* qmaxROC = 
201           fTPCdataQA->GetMaxCharge()->GetCalROC(iSec); 
202
203         const Int_t nRows = occupancyROC->GetNrows(); 
204         for (Int_t iRow = 0; iRow < nRows; iRow++) {
205
206           Int_t helpSector = iSec;
207           if(iRow>=64)
208             helpSector += 36; // OROC (long pads)
209
210           const Int_t nPads = occupancyROC->GetNPads(iRow); 
211           for (Int_t iPad = 0; iPad < nPads; iPad++) {
212         
213             fHistRawsOccupancy->Fill(occupancyROC->GetValue(iRow, iPad));
214             hNorm72->Fill(iSec);
215             fHistRawsOccupancyVsSector
216               ->Fill(iSec, occupancyROC->GetValue(iRow, iPad));
217
218             const Int_t nClusters = nclusterROC->GetValue(iRow, iPad);
219             
220             if(nClusters>0) {
221               
222               fHistRawsNClustersPerEventVsSector->Fill(iSec, nClusters);
223               hNorm108->Fill(helpSector, nClusters);
224               fHistRawsQVsSector->Fill(helpSector, 
225                                        nClusters*qROC->GetValue(iRow, iPad));
226               fHistRawsQmaxVsSector->Fill(helpSector, 
227                                           nClusters*qmaxROC->GetValue(iRow, iPad));
228             }
229           }
230         }
231       } // end loop over sectors
232       
233       // Normalize histograms
234       fHistRawsOccupancyVsSector->Divide(hNorm72);
235       fHistRawsNClustersPerEventVsSector->Scale(1.0/Float_t(eventCounter));
236       fHistRawsQVsSector->Divide(hNorm108);
237       fHistRawsQmaxVsSector->Divide(hNorm108);
238       delete hNorm72;
239       delete hNorm108;
240
241     }
242   }
243
244   AliQAChecker::Instance()->Run(AliQA::kTPC, task, list) ;  
245 }
246
247 //____________________________________________________________________________ 
248 void AliTPCQADataMakerRec::InitESDs()
249 {
250   //create ESDs histograms in ESDs subdir  
251   fHistESDclusters = 
252     new TH1F("hESDclusters", "N TPC clusters per track; N clusters; Counts",
253              160, 0, 160);
254   fHistESDclusters->Sumw2();
255   Add2ESDsList(fHistESDclusters, 0);
256
257   fHistESDratio = 
258     new TH1F("hESDratio", "Ratio: TPC clusters / findable; Ratio: cluster/findable; Counts",
259              100, 0, 1);
260   fHistESDratio->Sumw2();
261   Add2ESDsList(fHistESDratio, 1);
262   
263   fHistESDpt = 
264     new TH1F("hESDpt", "P_{T} distribution; p_{T} [GeV/c]; Counts",
265              50, 0, 5);
266   fHistESDpt->Sumw2();
267   Add2ESDsList(fHistESDpt, 2);
268 }
269
270 //____________________________________________________________________________ 
271 void AliTPCQADataMakerRec::InitRaws()
272 {
273   //
274   // Adding the raw 
275   //  
276
277   // Modified: 7/7 - 2008
278   // Laurent Aphecetche pointed out that the mapping was read from file
279   // for each event, so now we read in the map here and set if for 
280   // the raw data qa
281   fTPCdataQA = new AliTPCdataQA();
282   LoadMaps(); // Load Altro maps
283   fTPCdataQA->SetAltroMapping(fMapping); // set Altro mapping
284   fTPCdataQA->SetRangeTime(100, 920); // set time bin interval 
285   Add2RawsList(fTPCdataQA, 0); // This is used by the AMORE monitoring
286
287   fHistRawsOccupancy = 
288     new TH1F("hRawsOccupancy", "Occupancy (all pads); Occupancy; Counts",
289              100, 0, 1);
290   fHistRawsOccupancy->Sumw2();
291   Add2RawsList(fHistRawsOccupancy, 1);
292
293   fHistRawsOccupancyVsSector = 
294     new TH1F("hRawsOccupancyVsSector", "Occupancy vs sector; Sector; Occupancy",
295              72, 0, 72);
296   fHistRawsOccupancyVsSector->Sumw2();
297   Add2RawsList(fHistRawsOccupancyVsSector, 2);
298
299   fHistRawsNClustersPerEventVsSector = 
300     new TH1F("hRawsNClustersPerEventVsSector", "Nclusters per event vs sector; Sector; Nclusters per event",
301              72, 0, 72);
302   fHistRawsNClustersPerEventVsSector->Sumw2();
303   Add2RawsList(fHistRawsNClustersPerEventVsSector, 3);
304   
305   fHistRawsQVsSector = 
306     new TH1F("hRawsQVsSector", "<Q> vs sector (OROC med: 36-71, long: 72-107); Sector; <Q>",
307              108, 0, 108);
308   fHistRawsQVsSector->Sumw2();
309   Add2RawsList(fHistRawsQVsSector, 4);
310
311   fHistRawsQmaxVsSector = 
312     new TH1F("hRawsQmaxVsSector", "<Qmax> vs sector (OROC med: 36-71, long: 72-107); Sector; <Qmax>",
313              108, 0, 108);
314   fHistRawsQmaxVsSector->Sumw2();
315   Add2RawsList(fHistRawsQmaxVsSector, 5);
316 }
317
318 //____________________________________________________________________________ 
319 void AliTPCQADataMakerRec::InitRecPoints()
320 {
321   fHistRecPointsQmaxShort = 
322     new TH1F("hRecPointsQmaxShort", "Qmax distrbution (short pads); Qmax; Counts",
323              100, 0, 300);
324   fHistRecPointsQmaxShort->Sumw2();
325   Add2RecPointsList(fHistRecPointsQmaxShort, 0);
326
327   fHistRecPointsQmaxMedium = 
328     new TH1F("hRecPointsQmaxMedium", "Qmax distrbution (medium pads); Qmax; Counts",
329              100, 0, 300);
330   fHistRecPointsQmaxMedium->Sumw2();
331   Add2RecPointsList(fHistRecPointsQmaxMedium, 1);
332
333   fHistRecPointsQmaxLong = 
334     new TH1F("hRecPointsQmaxLong", "Qmax distrbution (long pads); Qmax; Counts",
335              100, 0, 300);
336   fHistRecPointsQmaxLong->Sumw2();
337   Add2RecPointsList(fHistRecPointsQmaxLong, 2);
338
339   fHistRecPointsQShort = 
340     new TH1F("hRecPointsQShort", "Q distrbution (short pads); Q; Counts",
341              100, 0, 2000);
342   fHistRecPointsQShort->Sumw2();
343   Add2RecPointsList(fHistRecPointsQShort, 3);
344
345   fHistRecPointsQMedium = 
346     new TH1F("hRecPointsQMedium", "Q distrbution (medium pads); Q; Counts",
347              100, 0, 2000);
348   fHistRecPointsQMedium->Sumw2();
349   Add2RecPointsList(fHistRecPointsQMedium, 4);
350
351   fHistRecPointsQLong = 
352     new TH1F("hRecPointsQLong", "Q distrbution (long pads); Q; Counts",
353              100, 0, 2000);
354   fHistRecPointsQLong->Sumw2();
355   Add2RecPointsList(fHistRecPointsQLong, 5);
356
357   fHistRecPointsRow = 
358     new TH1F("hRecPointsRow", "Clusters per row; Row; Counts",
359              159, 0, 159);
360   fHistRecPointsRow->Sumw2();
361   Add2RecPointsList(fHistRecPointsRow, 6);
362 }
363
364 //____________________________________________________________________________
365 void AliTPCQADataMakerRec::MakeESDs(AliESDEvent * esd)
366 {
367   // make QA data from ESDs
368   
369   const Int_t nESDTracks = esd->GetNumberOfTracks();
370   Int_t nTPCtracks = 0; 
371   for(Int_t i = 0; i < nESDTracks; i++) {
372     
373     AliESDtrack * track = esd->GetTrack(i);
374     
375     if ((track->GetStatus() & AliESDtrack::kTPCrefit)==0)
376       continue;
377     
378     nTPCtracks++;
379     
380     Int_t nTPCclusters         = track->GetTPCNcls();
381     Int_t nTPCclustersFindable = track->GetTPCNclsF();
382     if ( nTPCclustersFindable<=0) continue;
383     fHistESDclusters->Fill(nTPCclusters);
384     fHistESDratio->Fill(Float_t(nTPCclusters)/Float_t(nTPCclustersFindable));
385     fHistESDpt->Fill(track->Pt()); 
386   }
387 }
388
389 //____________________________________________________________________________
390 void AliTPCQADataMakerRec::MakeRaws(AliRawReader* rawReader)
391 {
392   //
393   // To make QA for the RAW data we use the TPC Calibration framework 
394   // to handle the data and then in the end extract the data
395   //
396         rawReader->Reset() ; 
397   fTPCdataQA->ProcessEvent(rawReader);  
398 }
399
400 //____________________________________________________________________________
401 void AliTPCQADataMakerRec::MakeRecPoints(TTree* recTree)
402 {
403   AliTPCClustersRow *clrow = new AliTPCClustersRow();
404   clrow->SetClass("AliTPCclusterMI");
405   clrow->SetArray(0);
406   clrow->GetArray()->ExpandCreateFast(10000);
407
408   TBranch* branch = recTree->GetBranch("Segment");
409   branch->SetAddress(&clrow);
410
411   const Int_t nEntries = Int_t(recTree->GetEntries());
412   for (Int_t i = 0; i < nEntries; i++) {
413     
414     branch->GetEntry(i);
415     
416     const Int_t nClusters = clrow->GetArray()->GetEntriesFast();
417     for (Int_t icl=0; icl < nClusters; icl++){
418       
419       AliTPCclusterMI* cluster = 
420         (AliTPCclusterMI*)clrow->GetArray()->At(icl);
421       
422       Float_t Qmax = cluster->GetMax();
423       Float_t Q    = cluster->GetQ();
424       Int_t   row  = cluster->GetRow();
425
426       if(cluster->GetDetector()<36) { // IROC (short pads)
427
428         fHistRecPointsQmaxShort->Fill(Qmax);
429         fHistRecPointsQShort->Fill(Q);
430       } else { // OROC (medium and long pads)
431         row += 63;
432         if(cluster->GetRow()<64) { // medium pads
433
434           fHistRecPointsQmaxMedium->Fill(Qmax);
435           fHistRecPointsQMedium->Fill(Q);
436         } else { // long pads
437
438           fHistRecPointsQmaxLong->Fill(Qmax);
439           fHistRecPointsQLong->Fill(Q);
440         }
441       }
442       
443       fHistRecPointsRow->Fill(row);
444     } // end loop over clusters
445   } // end loop over tree
446
447   delete clrow;
448 }
449
450 //____________________________________________________________________________
451 void AliTPCQADataMakerRec::LoadMaps()
452 {
453   TString path = gSystem->Getenv("ALICE_ROOT");
454   path += "/TPC/mapping/Patch";
455
456   for(Int_t i = 0; i < 6; i++) {
457     TString path2 = path;
458     path2 += i;
459     path2 += ".data";
460     fMapping[i] = new AliTPCAltroMapping(path2.Data());
461   }
462 }
463