1.The QA data created on demand according to the event species at filling time. 2...
[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   // Check id histograms already created for this Event Specie
382   if ( ! GetESDsData(KClusters) )
383     InitESDs() ;
384   
385   const Int_t nESDTracks = esd->GetNumberOfTracks();
386   Int_t nTPCtracks = 0; 
387   for(Int_t i = 0; i < nESDTracks; i++) {
388     
389     AliESDtrack * track = esd->GetTrack(i);
390     
391     if ((track->GetStatus() & AliESDtrack::kTPCrefit)==0)
392       continue;
393     
394     nTPCtracks++;
395     
396     Int_t nTPCclusters         = track->GetTPCNcls();
397     Int_t nTPCclustersFindable = track->GetTPCNclsF();
398     if ( nTPCclustersFindable<=0) continue;
399     GetESDsData(KClusters)->Fill(nTPCclusters);
400     GetESDsData(kRatio)->Fill(Float_t(nTPCclusters)/Float_t(nTPCclustersFindable));
401     GetESDsData(kPt)->Fill(track->Pt()); 
402   }
403 }
404
405 //____________________________________________________________________________
406 void AliTPCQADataMakerRec::MakeRaws(AliRawReader* rawReader)
407 {
408   //
409   // To make QA for the RAW data we use the TPC Calibration framework 
410   // to handle the data and then in the end extract the data
411   //
412   if ( ! GetRawsData(kOccupancy) )
413     InitRaws() ;
414   
415   rawReader->Reset() ; 
416   fTPCdataQA[AliRecoParam::AConvert(fEventSpecie)]->ProcessEvent(rawReader);  
417  }
418
419 //____________________________________________________________________________
420 void AliTPCQADataMakerRec::MakeDigits(TTree* digitTree)
421 {
422
423   // Check id histograms already created for this Event Specie
424   if ( ! GetDigitsData(kDigitsADC) )
425     InitDigits() ;
426   
427   TBranch* branch = digitTree->GetBranch("Segment");
428   AliSimDigits* digArray = 0;
429   branch->SetAddress(&digArray);
430   
431   Int_t nEntries = Int_t(digitTree->GetEntries());
432   
433   for (Int_t n = 0; n < nEntries; n++) {
434     
435     digitTree->GetEvent(n);
436     
437     if (digArray->First())
438       do {
439         Float_t dig = digArray->CurrentDigit();
440         
441         GetDigitsData(kDigitsADC)->Fill(dig);
442       } while (digArray->Next());    
443   }
444 }
445
446 //____________________________________________________________________________
447 void AliTPCQADataMakerRec::MakeRecPoints(TTree* recTree)
448 {
449   // Check id histograms already created for this Event Specie
450   if ( ! GetRecPointsData(kQmaxShort) )
451     InitRecPoints() ;
452
453   AliTPCClustersRow *clrow = new AliTPCClustersRow();
454   clrow->SetClass("AliTPCclusterMI");
455   clrow->SetArray(0);
456   clrow->GetArray()->ExpandCreateFast(10000);
457
458   TBranch* branch = recTree->GetBranch("Segment");
459   branch->SetAddress(&clrow);
460
461   const Int_t nEntries = Int_t(recTree->GetEntries());
462   for (Int_t i = 0; i < nEntries; i++) {
463     
464     branch->GetEntry(i);
465     
466     const Int_t nClusters = clrow->GetArray()->GetEntriesFast();
467     for (Int_t icl=0; icl < nClusters; icl++){
468       
469       AliTPCclusterMI* cluster = 
470         (AliTPCclusterMI*)clrow->GetArray()->At(icl);
471       
472       Float_t Qmax = cluster->GetMax();
473       Float_t Q    = cluster->GetQ();
474       Int_t   row  = cluster->GetRow();
475
476       if(cluster->GetDetector()<36) { // IROC (short pads)
477
478         GetRecPointsData(kQmaxShort)->Fill(Qmax);
479         GetRecPointsData(kQShort)->Fill(Q);
480       } else { // OROC (medium and long pads)
481         row += 63;
482         if(cluster->GetRow()<64) { // medium pads
483
484           GetRecPointsData(kQmaxMedium)->Fill(Qmax);
485           GetRecPointsData(kQMedium)->Fill(Q);
486         } else { // long pads
487
488           GetRecPointsData(kQmaxLong)->Fill(Qmax);
489           GetRecPointsData(kQLong)->Fill(Q);
490         }
491       }
492       
493       GetRecPointsData(kRow)->Fill(row);
494     } // end loop over clusters
495   } // end loop over tree
496
497   delete clrow;
498 }
499
500 //____________________________________________________________________________
501 void AliTPCQADataMakerRec::LoadMaps()
502 {
503   TString path = gSystem->Getenv("ALICE_ROOT");
504   path += "/TPC/mapping/Patch";
505
506   for(Int_t i = 0; i < 6; i++) {
507     TString path2 = path;
508     path2 += i;
509     path2 += ".data";
510     fMapping[i] = new AliTPCAltroMapping(path2.Data());
511   }
512 }
513