Added QA for digits during reconstruction (Yves)
[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 #include "AliSimDigits.h"
67
68 ClassImp(AliTPCQADataMakerRec)
69
70 //____________________________________________________________________________ 
71 AliTPCQADataMakerRec::AliTPCQADataMakerRec() : 
72   AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kTPC), 
73                     "TPC Rec Quality Assurance Data Maker"),
74   fTPCdataQA(NULL)
75 {
76   // ctor
77   fTPCdataQA = new AliTPCdataQA*[AliRecoParam::kNSpecies] ;
78   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
79     fTPCdataQA[specie] = NULL ; 
80   
81   for(Int_t i = 0; i < 6; i++)
82     fMapping[i] = 0;
83 }
84
85 //____________________________________________________________________________ 
86 AliTPCQADataMakerRec::AliTPCQADataMakerRec(const AliTPCQADataMakerRec& qadm) :
87   AliQADataMakerRec(),
88   fTPCdataQA(NULL)
89 {
90   //copy ctor 
91   // Does not copy the calibration object, instead InitRaws have to be
92   // called again
93   SetName((const char*)qadm.GetName()) ; 
94   SetTitle((const char*)qadm.GetTitle()); 
95
96   fTPCdataQA = new AliTPCdataQA*[AliRecoParam::kNSpecies] ;
97   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
98     fTPCdataQA[specie] = NULL ; 
99   
100   for(Int_t i = 0; i < 6; i++)
101     fMapping[i] = 0;
102
103   //
104   // Associate class histogram objects to the copies in the list
105   // Could also be done with the indexes
106   //
107
108 }
109
110 //__________________________________________________________________
111 AliTPCQADataMakerRec& AliTPCQADataMakerRec::operator = (const AliTPCQADataMakerRec& qadm )
112 {
113   // Equal operator.
114   this->~AliTPCQADataMakerRec();
115   new(this) AliTPCQADataMakerRec(qadm);
116   return *this;
117 }
118
119 //__________________________________________________________________
120 AliTPCQADataMakerRec::~AliTPCQADataMakerRec()
121 {
122   // Destructor
123   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) 
124     if ( fTPCdataQA[specie] != NULL )
125       delete fTPCdataQA[specie] ; 
126   delete[] fTPCdataQA; 
127
128   for(Int_t i = 0; i < 6; i++) 
129     delete fMapping[i];
130 }
131  
132 //____________________________________________________________________________ 
133 void AliTPCQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
134 {
135   //Detector specific actions at end of cycle
136   
137   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
138     if ( !AliQAv1::Instance()->IsEventSpecieSet(specie) ) 
139       continue ; 
140     if(fTPCdataQA[specie] != NULL) { // do the final step of the QA for Raw data
141
142       fTPCdataQA[specie]->Analyse(); // 31/1-08 Analyse is now protected against
143                            //         RAW data files with no TPC data
144
145       SetEventSpecie(AliRecoParam::ConvertIndex(specie)) ; 
146       TH1F * histRawsOccupancy                 = (TH1F*)GetRawsData(kOccupancy) ;
147       TH1F * histRawsOccupancyVsSector         = (TH1F*)GetRawsData(kOccupancyVsSector) ;
148       TH1F * histRawsNClustersPerEventVsSector = (TH1F*)GetRawsData(kNClustersPerEventVsSector) ;
149       TH1F * histRawsQVsSector                 = (TH1F*)GetRawsData(kQVsSector) ;
150       TH1F * histRawsQmaxVsSector              = (TH1F*)GetRawsData(kQmaxVsSector) ;
151       if ( !histRawsOccupancy ||
152           !histRawsOccupancyVsSector ||
153           !histRawsNClustersPerEventVsSector ||
154           !histRawsQVsSector ||
155           !histRawsQmaxVsSector) {
156         AliError("Something very wrong here, corrupted memory ?????. Please check\n") ; 
157         continue ; 
158       }
159         
160       //Add2RawsList(fTPCdataQA, 0);
161       // get the histograms and add them to the output
162       // 31/8-08 Histogram is only added if the Calibration class 
163       //         receives TPC data 
164       const Int_t eventCounter = fTPCdataQA[specie]->GetEventCounter();
165       if(eventCounter>0) { // some TPC data has been processed
166
167         // Reset histograms and refill them 
168         histRawsOccupancy->Reset();
169         histRawsOccupancyVsSector->Reset();
170         histRawsNClustersPerEventVsSector->Reset();
171         histRawsQVsSector->Reset();
172         histRawsQmaxVsSector->Reset();
173       
174         TH1F* hNorm72 = new TH1F("hNorm72", "histogram to normalize 72 sectors",
175                                  72, 0, 72);
176         hNorm72->Sumw2();
177         TH1F* hNorm108 = new TH1F("hNorm108", "histogram to normalize 108 sectors (medium and long pads are split up)",
178                                   108, 0, 108);
179         hNorm108->Sumw2();
180
181         for (Int_t iSec = 0; iSec < 72; iSec++) {
182         
183           AliTPCCalROC* occupancyROC = 
184           fTPCdataQA[specie]->GetNoThreshold()->GetCalROC(iSec); 
185           AliTPCCalROC* nclusterROC = 
186           fTPCdataQA[specie]->GetNLocalMaxima()->GetCalROC(iSec); 
187           AliTPCCalROC* qROC = 
188           fTPCdataQA[specie]->GetMeanCharge()->GetCalROC(iSec); 
189           AliTPCCalROC* qmaxROC = 
190           fTPCdataQA[specie]->GetMaxCharge()->GetCalROC(iSec); 
191
192           const Int_t nRows = occupancyROC->GetNrows(); 
193           for (Int_t iRow = 0; iRow < nRows; iRow++) {
194
195             Int_t helpSector = iSec;
196             if(iRow>=64)
197               helpSector += 36; // OROC (long pads)
198
199             const Int_t nPads = occupancyROC->GetNPads(iRow); 
200             for (Int_t iPad = 0; iPad < nPads; iPad++) {
201         
202               histRawsOccupancy->Fill(occupancyROC->GetValue(iRow, iPad));
203               hNorm72->Fill(iSec);
204               histRawsOccupancyVsSector
205               ->Fill(iSec, occupancyROC->GetValue(iRow, iPad));
206
207               const Int_t nClusters = TMath::Nint(nclusterROC->GetValue(iRow, iPad));
208             
209               if(nClusters>0) {
210               
211                 histRawsNClustersPerEventVsSector->Fill(iSec, nClusters);
212                 hNorm108->Fill(helpSector, nClusters);
213                 histRawsQVsSector->Fill(helpSector, 
214                                        nClusters*qROC->GetValue(iRow, iPad));
215                 histRawsQmaxVsSector->Fill(helpSector, 
216                                           nClusters*qmaxROC->GetValue(iRow, iPad));
217               }
218             }
219           }
220         } // end loop over sectors
221       
222         // Normalize histograms
223         histRawsOccupancyVsSector->Divide(hNorm72);
224         histRawsNClustersPerEventVsSector->Scale(1.0/Float_t(eventCounter));
225         histRawsQVsSector->Divide(hNorm108);
226         histRawsQmaxVsSector->Divide(hNorm108);
227         delete hNorm72;
228         delete hNorm108;
229       }
230     }
231   }
232   AliQAChecker::Instance()->Run(AliQAv1::kTPC, task, list) ;  
233 }
234
235 //____________________________________________________________________________ 
236 void AliTPCQADataMakerRec::InitESDs()
237 {
238   //create ESDs histograms in ESDs subdir  
239   const Bool_t expert   = kTRUE ; 
240   const Bool_t image    = kTRUE ; 
241   
242   TH1F * histESDclusters = 
243     new TH1F("hESDclusters", "N TPC clusters per track; N clusters; Counts",
244              160, 0, 160);
245   histESDclusters->Sumw2();
246   Add2ESDsList(histESDclusters, KClusters, !expert, image);
247
248   TH1F * histESDratio = 
249     new TH1F("hESDratio", "Ratio: TPC clusters / findable; Ratio: cluster/findable; Counts",
250              100, 0, 1);
251   histESDratio->Sumw2();
252   Add2ESDsList(histESDratio, kRatio, !expert, image);
253   
254   TH1F * histESDpt = 
255     new TH1F("hESDpt", "P_{T} distribution; p_{T} [GeV/c]; Counts",
256              50, 0, 5);
257   histESDpt->Sumw2();
258   Add2ESDsList(histESDpt, kPt, !expert, image);
259 }
260
261 //____________________________________________________________________________ 
262 void AliTPCQADataMakerRec::InitRaws()
263 {
264   //
265   // Adding the raw 
266   //  
267
268   // Modified: 7/7 - 2008
269   // Laurent Aphecetche pointed out that the mapping was read from file
270   // for each event, so now we read in the map here and set if for 
271   // the raw data qa
272   const Bool_t expert   = kTRUE ; 
273   const Bool_t saveCorr = kTRUE ; 
274   const Bool_t image    = kTRUE ; 
275   
276   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
277     fTPCdataQA[specie] = new AliTPCdataQA(AliRecoParam::ConvertIndex(specie));
278     LoadMaps(); // Load Altro maps
279     fTPCdataQA[specie]->SetAltroMapping(fMapping); // set Altro mapping
280     fTPCdataQA[specie]->SetRangeTime(100, 920); // set time bin interval 
281 //    Add2RawsList(fTPCdataQA, kTPCdataQ, !expert, image, !saveCorrA); // This is used by the AMORE monitoring <------- THIS WILL FAIL (YS)
282   }
283
284   TH1F * histRawsOccupancy = 
285     new TH1F("hRawsOccupancy", "Occupancy (all pads); Occupancy; Counts",
286              100, 0, 1);
287   histRawsOccupancy->Sumw2();
288   Add2RawsList(histRawsOccupancy, kOccupancy, !expert, image, !saveCorr);
289   
290   TH1F * histRawsOccupancyVsSector = 
291     new TH1F("hRawsOccupancyVsSector", "Occupancy vs sector; Sector; Occupancy",
292              72, 0, 72);
293   histRawsOccupancyVsSector->Sumw2();
294   Add2RawsList(histRawsOccupancyVsSector, kOccupancyVsSector, !expert, image, !saveCorr);
295
296   TH1F * histRawsNClustersPerEventVsSector = 
297     new TH1F("hRawsNClustersPerEventVsSector", "Nclusters per event vs sector; Sector; Nclusters per event",
298              72, 0, 72);
299   histRawsNClustersPerEventVsSector->Sumw2();
300   Add2RawsList(histRawsNClustersPerEventVsSector, kNClustersPerEventVsSector, !expert, image, !saveCorr);
301   
302   TH1F * histRawsQVsSector = 
303     new TH1F("hRawsQVsSector", "<Q> vs sector (OROC med: 36-71, long: 72-107); Sector; <Q>",
304              108, 0, 108);
305   histRawsQVsSector->Sumw2();
306   Add2RawsList(histRawsQVsSector, kQVsSector, !expert, image, !saveCorr);
307
308   TH1F * histRawsQmaxVsSector = 
309     new TH1F("hRawsQmaxVsSector", "<Qmax> vs sector (OROC med: 36-71, long: 72-107); Sector; <Qmax>",
310              108, 0, 108);
311   histRawsQmaxVsSector->Sumw2();
312   Add2RawsList(histRawsQmaxVsSector, kQmaxVsSector, !expert, image, !saveCorr);
313 }
314
315 //____________________________________________________________________________ 
316 void AliTPCQADataMakerRec::InitDigits()
317 {
318   const Bool_t expert   = kTRUE ; 
319   const Bool_t image    = kTRUE ; 
320   TH1F * histDigitsADC = 
321     new TH1F("hDigitsADC", "Digit ADC distribution; ADC; Counts",
322              1000, 0, 1000);
323   histDigitsADC->Sumw2();
324   Add2DigitsList(histDigitsADC, kDigitsADC, !expert, image);
325 }
326
327 //____________________________________________________________________________ 
328 void AliTPCQADataMakerRec::InitRecPoints()
329 {
330   const Bool_t expert   = kTRUE ; 
331   const Bool_t image    = kTRUE ; 
332   
333   TH1F * histRecPointsQmaxShort = 
334     new TH1F("hRecPointsQmaxShort", "Qmax distrbution (short pads); Qmax; Counts",
335              100, 0, 300);
336   histRecPointsQmaxShort->Sumw2();
337   Add2RecPointsList(histRecPointsQmaxShort, kQmaxShort, !expert, image);
338
339   TH1F * histRecPointsQmaxMedium = 
340     new TH1F("hRecPointsQmaxMedium", "Qmax distrbution (medium pads); Qmax; Counts",
341              100, 0, 300);
342   histRecPointsQmaxMedium->Sumw2();
343   Add2RecPointsList(histRecPointsQmaxMedium, kQmaxMedium, !expert, image);
344
345   TH1F * histRecPointsQmaxLong = 
346     new TH1F("hRecPointsQmaxLong", "Qmax distrbution (long pads); Qmax; Counts",
347              100, 0, 300);
348   histRecPointsQmaxLong->Sumw2();
349   Add2RecPointsList(histRecPointsQmaxLong, kQmaxLong, !expert, image);
350
351   TH1F * histRecPointsQShort = 
352     new TH1F("hRecPointsQShort", "Q distrbution (short pads); Q; Counts",
353              100, 0, 2000);
354   histRecPointsQShort->Sumw2();
355   Add2RecPointsList(histRecPointsQShort, kQShort, !expert, image);
356
357   TH1F * histRecPointsQMedium = 
358     new TH1F("hRecPointsQMedium", "Q distrbution (medium pads); Q; Counts",
359              100, 0, 2000);
360   histRecPointsQMedium->Sumw2();
361   Add2RecPointsList(histRecPointsQMedium, kQMedium, !expert, image);
362
363   TH1F * histRecPointsQLong = 
364     new TH1F("hRecPointsQLong", "Q distrbution (long pads); Q; Counts",
365              100, 0, 2000);
366   histRecPointsQLong->Sumw2();
367   Add2RecPointsList(histRecPointsQLong, kQLong, !expert, image);
368
369   TH1F * histRecPointsRow = 
370     new TH1F("hRecPointsRow", "Clusters per row; Row; Counts",
371              159, 0, 159);
372   histRecPointsRow->Sumw2();
373   Add2RecPointsList(histRecPointsRow, kRow, !expert, image);
374 }
375
376 //____________________________________________________________________________
377 void AliTPCQADataMakerRec::MakeESDs(AliESDEvent * esd)
378 {
379   // make QA data from ESDs
380   
381   const Int_t nESDTracks = esd->GetNumberOfTracks();
382   Int_t nTPCtracks = 0; 
383   for(Int_t i = 0; i < nESDTracks; i++) {
384     
385     AliESDtrack * track = esd->GetTrack(i);
386     
387     if ((track->GetStatus() & AliESDtrack::kTPCrefit)==0)
388       continue;
389     
390     nTPCtracks++;
391     
392     Int_t nTPCclusters         = track->GetTPCNcls();
393     Int_t nTPCclustersFindable = track->GetTPCNclsF();
394     if ( nTPCclustersFindable<=0) continue;
395     GetESDsData(KClusters)->Fill(nTPCclusters);
396     GetESDsData(kRatio)->Fill(Float_t(nTPCclusters)/Float_t(nTPCclustersFindable));
397     GetESDsData(kPt)->Fill(track->Pt()); 
398   }
399 }
400
401 //____________________________________________________________________________
402 void AliTPCQADataMakerRec::MakeRaws(AliRawReader* rawReader)
403 {
404   //
405   // To make QA for the RAW data we use the TPC Calibration framework 
406   // to handle the data and then in the end extract the data
407   //
408         rawReader->Reset() ; 
409   fTPCdataQA[AliRecoParam::AConvert(fEventSpecie)]->ProcessEvent(rawReader);  
410  }
411
412 //____________________________________________________________________________
413 void AliTPCQADataMakerRec::MakeDigits(TTree* digitTree)
414 {
415   TBranch* branch = digitTree->GetBranch("Segment");
416   AliSimDigits* digArray = 0;
417   branch->SetAddress(&digArray);
418   
419   Int_t nEntries = Int_t(digitTree->GetEntries());
420   
421   for (Int_t n = 0; n < nEntries; n++) {
422     
423     digitTree->GetEvent(n);
424     
425     if (digArray->First())
426       do {
427         Float_t dig = digArray->CurrentDigit();
428         
429         GetDigitsData(kDigitsADC)->Fill(dig);
430       } while (digArray->Next());    
431   }
432 }
433
434 //____________________________________________________________________________
435 void AliTPCQADataMakerRec::MakeRecPoints(TTree* recTree)
436 {
437   AliTPCClustersRow *clrow = new AliTPCClustersRow();
438   clrow->SetClass("AliTPCclusterMI");
439   clrow->SetArray(0);
440   clrow->GetArray()->ExpandCreateFast(10000);
441
442   TBranch* branch = recTree->GetBranch("Segment");
443   branch->SetAddress(&clrow);
444
445   const Int_t nEntries = Int_t(recTree->GetEntries());
446   for (Int_t i = 0; i < nEntries; i++) {
447     
448     branch->GetEntry(i);
449     
450     const Int_t nClusters = clrow->GetArray()->GetEntriesFast();
451     for (Int_t icl=0; icl < nClusters; icl++){
452       
453       AliTPCclusterMI* cluster = 
454         (AliTPCclusterMI*)clrow->GetArray()->At(icl);
455       
456       Float_t Qmax = cluster->GetMax();
457       Float_t Q    = cluster->GetQ();
458       Int_t   row  = cluster->GetRow();
459
460       if(cluster->GetDetector()<36) { // IROC (short pads)
461
462         GetRecPointsData(kQmaxShort)->Fill(Qmax);
463         GetRecPointsData(kQShort)->Fill(Q);
464       } else { // OROC (medium and long pads)
465         row += 63;
466         if(cluster->GetRow()<64) { // medium pads
467
468           GetRecPointsData(kQmaxMedium)->Fill(Qmax);
469           GetRecPointsData(kQMedium)->Fill(Q);
470         } else { // long pads
471
472           GetRecPointsData(kQmaxLong)->Fill(Qmax);
473           GetRecPointsData(kQLong)->Fill(Q);
474         }
475       }
476       
477       GetRecPointsData(kRow)->Fill(row);
478     } // end loop over clusters
479   } // end loop over tree
480
481   delete clrow;
482 }
483
484 //____________________________________________________________________________
485 void AliTPCQADataMakerRec::LoadMaps()
486 {
487   TString path = gSystem->Getenv("ALICE_ROOT");
488   path += "/TPC/mapping/Patch";
489
490   for(Int_t i = 0; i < 6; i++) {
491     TString path2 = path;
492     path2 += i;
493     path2 += ".data";
494     fMapping[i] = new AliTPCAltroMapping(path2.Data());
495   }
496 }
497