1.The QA data created on demand according to the event species at filling time. 2...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALQADataMakerRec.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 Based on the QA code for PHOS written by Yves Schutz July 2007
17
18 Authors:  J.Klay (Cal Poly) May 2008
19           S. Salur LBL April 2008
20
21 */
22
23 // --- ROOT system ---
24 #include <TClonesArray.h>
25 #include <TFile.h> 
26 #include <TH1F.h> 
27 #include <TH1I.h> 
28 #include <TH2F.h> 
29 #include <TProfile.h> 
30
31 // --- Standard library ---
32
33
34 // --- AliRoot header files ---
35 #include "AliESDCaloCluster.h"
36 #include "AliESDCaloCells.h"
37 #include "AliESDEvent.h"
38 #include "AliLog.h"
39 #include "AliEMCALQADataMakerRec.h"
40 #include "AliQAChecker.h"
41 #include "AliEMCALDigit.h" 
42 #include "AliEMCALRecPoint.h" 
43 #include "AliEMCALRawUtils.h"
44 #include "AliEMCALReconstructor.h"
45 #include "AliEMCALRecParam.h"
46 #include "AliRawReader.h"
47 #include "AliCaloRawStream.h"
48
49 ClassImp(AliEMCALQADataMakerRec)
50            
51 //____________________________________________________________________________ 
52   AliEMCALQADataMakerRec::AliEMCALQADataMakerRec() : 
53     AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kEMCAL), "EMCAL Quality Assurance Data Maker"),
54     fSuperModules(4) // FIXME!!! number of SuperModules; 4 for 2009; update default to 12 for later runs..
55 {
56   // ctor
57 }
58
59 //____________________________________________________________________________ 
60 AliEMCALQADataMakerRec::AliEMCALQADataMakerRec(const AliEMCALQADataMakerRec& qadm) :
61   AliQADataMakerRec(), fSuperModules()
62 {
63   //copy ctor 
64   SetName((const char*)qadm.GetName()) ; 
65   SetTitle((const char*)qadm.GetTitle()); 
66   fSuperModules = qadm.GetSuperModules();
67 }
68
69 //__________________________________________________________________
70 AliEMCALQADataMakerRec& AliEMCALQADataMakerRec::operator = (const AliEMCALQADataMakerRec& qadm )
71 {
72   // Equal operator.
73   this->~AliEMCALQADataMakerRec();
74   new(this) AliEMCALQADataMakerRec(qadm);
75   return *this;
76 }
77  
78 //____________________________________________________________________________ 
79 void AliEMCALQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
80 {
81   //Detector specific actions at end of cycle
82   // do the QA checking
83   AliQAChecker::Instance()->Run(AliQAv1::kEMCAL, task, list) ;  
84 }
85
86 //____________________________________________________________________________ 
87 void AliEMCALQADataMakerRec::InitESDs()
88 {
89   //Create histograms to controll ESD
90   const Bool_t expert   = kTRUE ; 
91   const Bool_t image    = kTRUE ; 
92   
93   TH1F * h1 = new TH1F("hESDCaloClusterE",  "ESDs CaloCluster energy in EMCAL;Energy [MeV];Counts",    200, 0., 20.) ; 
94   h1->Sumw2() ;
95   Add2ESDsList(h1, kESDCaloClusE, !expert, image)  ;                                                     
96
97   TH1I * h2 = new TH1I("hESDCaloClusterM", "ESDs CaloCluster multiplicity in EMCAL;# of Clusters;Entries", 100, 0,  100) ; 
98   h2->Sumw2() ;
99   Add2ESDsList(h2, kESDCaloClusM, !expert, image)  ;
100
101   TH1F * h3 = new TH1F("hESDCaloCellA",  "ESDs CaloCell amplitude in EMCAL;Energy [MeV];Counts",    500, 0., 250.) ; 
102   h3->Sumw2() ;
103   Add2ESDsList(h3, kESDCaloCellA, !expert, image)  ;  
104  
105   TH1I * h4 = new TH1I("hESDCaloCellM", "ESDs CaloCell multiplicity in EMCAL;# of Clusters;Entries", 200, 0,  1000) ; 
106   h4->Sumw2() ;
107   Add2ESDsList(h4, kESDCaloCellM, !expert, image) ;
108         
109 }
110
111 //____________________________________________________________________________ 
112 void AliEMCALQADataMakerRec::InitDigits()
113 {
114   // create Digits histograms in Digits subdir
115   const Bool_t expert   = kTRUE ; 
116   const Bool_t image    = kTRUE ; 
117   
118   TH1I * h0 = new TH1I("hEmcalDigits",    "Digits amplitude distribution in EMCAL;Amplitude [ADC counts];Counts",    500, 0, 500) ; 
119   h0->Sumw2() ;
120   Add2DigitsList(h0, 0, !expert, image) ;
121   TH1I * h1 = new TH1I("hEmcalDigitsMul", "Digits multiplicity distribution in EMCAL;# of Digits;Entries", 200, 0, 2000) ; 
122   h1->Sumw2() ;
123   Add2DigitsList(h1, 1, !expert, image) ;
124 }
125
126 //____________________________________________________________________________ 
127 void AliEMCALQADataMakerRec::InitRecPoints()
128 {
129   // create Reconstructed Points histograms in RecPoints subdir
130   const Bool_t expert   = kTRUE ; 
131   const Bool_t image    = kTRUE ; 
132   
133   TH1F* h0 = new TH1F("hEMCALRpE","EMCAL RecPoint energies;Energy [MeV];Counts",200, 0.,20.); //GeV
134   h0->Sumw2();
135   Add2RecPointsList(h0,kRecPE, !expert, image);
136
137   TH1I* h1 = new TH1I("hEMCALRpM","EMCAL RecPoint multiplicities;# of Clusters;Entries",100,0,100);
138   h1->Sumw2();
139   Add2RecPointsList(h1,kRecPM, !expert, image);
140
141   TH1I* h2 = new TH1I("hEMCALRpDigM","EMCAL RecPoint Digit Multiplicities;# of Digits;Entries",20,0,20);
142   h2->Sumw2();
143   Add2RecPointsList(h2,kRecPDigM, !expert, image);
144
145 }
146
147 //____________________________________________________________________________ 
148 void AliEMCALQADataMakerRec::InitRaws()
149 {
150   // create Raws histograms in Raws subdir
151    const Bool_t expert   = kTRUE ; 
152    const Bool_t saveCorr = kTRUE ; 
153    const Bool_t image    = kTRUE ; 
154
155   int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
156   int nTot = fSuperModules * nTowersPerSM; // max number of towers in all SuperModules
157
158   // counter info: number of channels per event (bins are SM index)
159   TProfile * h0 = new TProfile("hLowEmcalSupermodules", "Low Gain EMC: # of towers vs SuperMod;SM Id;# of towers",
160                                fSuperModules, -0.5, fSuperModules-0.5) ;
161   Add2RawsList(h0, kNsmodLG, !expert, image, !saveCorr) ;
162   TProfile * h1 = new TProfile("hHighEmcalSupermodules", "High Gain EMC: # of towers vs SuperMod;SM Id;# of towers",  
163                                fSuperModules, -0.5, fSuperModules-0.5) ; 
164   Add2RawsList(h1, kNsmodHG, !expert, image, !saveCorr) ;
165
166   // where did max sample occur? (bins are towers)
167   TProfile * h2 = new TProfile("hLowEmcalRawtime", "Low Gain EMC: Time at Max vs towerId;Tower Id;Time [ns]", 
168                                nTot, -0.5, nTot-0.5) ;
169   Add2RawsList(h2, kTimeLG, !expert, image, !saveCorr) ;
170   TProfile * h3 = new TProfile("hHighEmcalRawtime", "High Gain EMC: Time at Max vs towerId;Tower Id;Time [ns]", 
171                                nTot, -0.5, nTot-0.5) ;
172   Add2RawsList(h3, kTimeHG, !expert, image, !saveCorr) ;
173
174   // how much above pedestal was the max sample?  (bins are towers)
175   TProfile * h4 = new TProfile("hLowEmcalRawMaxMinusMin", "Low Gain EMC: Max - Min vs towerId;Tower Id;??", 
176                                nTot, -0.5, nTot-0.5) ;
177   Add2RawsList(h4, kSigLG, !expert, image, !saveCorr) ;
178   TProfile * h5 = new TProfile("hHighEmcalRawMaxMinusMin", "High Gain EMC: Max - Min vs towerId;Tower Id;??",
179                                nTot, -0.5, nTot-0.5) ;
180   Add2RawsList(h5, kSigHG, !expert, image, !saveCorr) ;
181
182   // total counter: channels per event
183   TH1I * h6 = new TH1I("hLowNtot", "Low Gain EMC: Total Number of found towers;# of Towers;Counts", 200, 0, nTot) ;
184   h6->Sumw2() ;
185   Add2RawsList(h6, kNtotLG, !expert, image, !saveCorr) ;
186   TH1I * h7 = new TH1I("hHighNtot", "High Gain EMC: Total Number of found towers;# of Towers;Counts", 200,0, nTot) ;
187   h7->Sumw2() ;
188   Add2RawsList(h7, kNtotHG, !expert, image, !saveCorr) ;
189
190   // pedestal (bins are towers)
191   TProfile * h8 = new TProfile("hLowEmcalRawPed", "Low Gain EMC: Pedestal vs towerId;Tower Id;Pedestal [ADC counts]", 
192                                nTot, -0.5, nTot-0.5) ;
193   Add2RawsList(h8, kPedLG, !expert, image, !saveCorr) ;
194   TProfile * h9 = new TProfile("hHighEmcalRawPed", "High Gain EMC: Pedestal vs towerId;Tower Id;Pedestal [ADC counts]",
195                                nTot, -0.5, nTot-0.5) ;
196   Add2RawsList(h9, kPedHG, !expert, image, !saveCorr) ;
197
198   // pedestal rms (standard dev = sqrt of variance estimator for pedestal) (bins are towers)
199   TProfile * h10 = new TProfile("hLowEmcalRawPedRMS", "Low Gain EMC: Pedestal RMS vs towerId;Tower Id;Width [ADC counts]", 
200                                 nTot, -0.5, nTot-0.5) ;
201   Add2RawsList(h10, kPedRMSLG, !expert, image, !saveCorr) ;
202   TProfile * h11 = new TProfile("hHighEmcalRawPedRMS", "High Gain EMC: Pedestal RMS vs towerId;Tower Id;Width [ADC counts]",
203                                 nTot, -0.5, nTot-0.5) ;
204   Add2RawsList(h11, kPedRMSHG, !expert, image, !saveCorr) ;
205   
206 }
207
208 //____________________________________________________________________________
209 void AliEMCALQADataMakerRec::MakeESDs(AliESDEvent * esd)
210 {
211   // make QA data from ESDs
212
213   
214   // Check id histograms already created for this Event Specie
215   if ( ! GetESDsData(kESDCaloClusE) )
216     InitESDs() ;
217   
218   Int_t nTot = 0 ; 
219   for ( Int_t index = 0; index < esd->GetNumberOfCaloClusters() ; index++ ) {
220     AliESDCaloCluster * clu = esd->GetCaloCluster(index) ;
221     if( clu->IsEMCAL() ) {
222       GetESDsData(kESDCaloClusE)->Fill(clu->E()) ;
223       nTot++ ;
224     } 
225   }
226   GetESDsData(kESDCaloClusM)->Fill(nTot) ;
227
228   //fill calo cells
229   AliESDCaloCells* cells = esd->GetEMCALCells();
230   GetESDsData(kESDCaloCellM)->Fill(cells->GetNumberOfCells()) ;
231
232   for ( Int_t index = 0; index < cells->GetNumberOfCells() ; index++ ) {
233     GetESDsData(kESDCaloCellA)->Fill(cells->GetAmplitude(index)) ;
234   }
235
236 }
237
238 //____________________________________________________________________________
239 void AliEMCALQADataMakerRec::MakeRaws(AliRawReader* rawReader)
240 {
241   //Fill prepared histograms with Raw digit properties
242
243   //Raw histogram filling not yet implemented
244   //
245   //Need to figure out how to get the info we want without having to
246   //actually run Raw2Digits twice.
247   //I suspect what we actually want is a raw digits method, not a true
248   //emcal raw data method, but this doesn't seem to be allowed in
249   //AliQADataMakerRec.h
250
251   // For now, to avoid redoing the expensive signal fits we just
252   // look at max vs min of the signal spextra, a la online usage in
253   // AliCaloCalibPedestal
254
255   // Check id histograms already created for this Event Specie
256   if ( ! GetRawsData(kSigLG) )
257     InitRaws() ;
258
259   rawReader->Reset() ;
260   AliCaloRawStream in(rawReader,"EMCAL"); 
261
262   // setup
263   int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
264   int nRows = 24; // number of rows per SuperModule
265   int sampleMin = 0; 
266   int sampleMax = 0x3ff; // 1023 = 10-bit range
267
268   // for the pedestal calculation
269   Bool_t selectPedestalSamples = kTRUE;
270   int firstPedestalSample = 0;
271   int lastPedestalSample = 15;
272
273   // SM counters; decl. should be safe, assuming we don't get more than 12 SuperModules..
274   int nTotalSMLG[12] = {0};
275   int nTotalSMHG[12] = {0};
276
277   // indices for the reading
278   int iSM = 0;
279   int sample = 0;
280   int gain = 0;
281   int time = 0;
282   // counters, on sample level
283   int i = 0; // the sample number in current event.
284   int max = sampleMin, min = sampleMax;//Use these for picking the pedestal
285   int maxTime = 0;
286
287   // for the pedestal calculation
288   int sampleSum = 0; // sum of samples
289   int squaredSampleSum = 0; // sum of samples squared
290   int nSum = 0; // number of samples in sum
291   // calc. quantities
292   double meanPed = 0, squaredMean = 0, rmsPed = 0;
293   
294   while (in.Next()) { // loop over input stream
295     sample = in.GetSignal(); //Get the adc signal
296     time = in.GetTime();
297     if (sample < min) { min = sample; }
298     if (sample > max) { 
299       max = sample;
300       maxTime = time;
301     }
302     i++;
303
304     // should we add it for the pedestal calculation?
305     if ( (firstPedestalSample<=time && time<=lastPedestalSample) || // sample time in range
306          !selectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all 
307       sampleSum += sample;
308       squaredSampleSum += sample*sample;
309       nSum++;
310     }
311
312     if ( i >= in.GetTimeLength()) {
313       // calculate pedesstal estimate: mean of possibly selected samples
314       if (nSum > 0) {
315         meanPed = sampleSum / (1.0 * nSum);
316         squaredMean = squaredSampleSum / (1.0 * nSum);
317         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
318         rmsPed = sqrt(squaredMean - meanPed*meanPed); 
319       }
320       else {
321         meanPed = 0;
322         squaredMean = 0;
323         rmsPed  = 0;
324       }
325
326       //If we're here then we're done with this tower
327       gain = -1; // init to not valid value
328       if ( in.IsLowGain() ) {
329         gain = 0;
330       }
331       else if ( in.IsHighGain() ) {
332         gain = 1;
333       }
334       
335       iSM = in.GetModule(); //The modules are numbered starting from 0
336
337       if (iSM>=0 && iSM<fSuperModules) { // valid module reading, can go on with filling
338
339         int towerId = iSM*nTowersPerSM + in.GetColumn()*nRows + in.GetRow();
340
341         if (gain == 0) { 
342           //fill the low gain histograms, and counters
343           nTotalSMLG[iSM]++; // one more channel found
344           GetRawsData(kSigLG)->Fill(towerId, max - min);
345           GetRawsData(kTimeLG)->Fill(towerId, maxTime);
346           if (nSum>0) { // only fill pedestal info in case it could be calculated
347             GetRawsData(kPedLG)->Fill(towerId, meanPed);
348             GetRawsData(kPedRMSLG)->Fill(towerId, rmsPed);
349           }
350         } // gain==0
351         else if (gain == 1) {           
352           //fill the high gain ones
353           nTotalSMHG[iSM]++; // one more channel found
354           GetRawsData(kSigHG)->Fill(towerId, max - min);
355           GetRawsData(kTimeHG)->Fill(towerId, maxTime);
356           if (nSum>0) { // only fill pedestal info in case it could be calculated
357             GetRawsData(kPedHG)->Fill(towerId, meanPed);
358             GetRawsData(kPedRMSHG)->Fill(towerId, rmsPed);
359           }
360         }
361       } // SM index OK
362
363       // reset counters
364       max = sampleMin; min = sampleMax;
365       maxTime = 0;
366       i = 0;
367       // also pedestal calc counters
368       sampleSum = 0; // sum of samples
369       squaredSampleSum = 0; // sum of samples squared
370       nSum = 0; // number of samples in sum
371     
372     }//End if, of channel
373    
374   }//end while, of stream
375
376   // let's also fill the SM and event counter histograms
377   int nTotalHG = 0;
378   int nTotalLG = 0;
379   for (iSM=0; iSM<fSuperModules; iSM++) {  
380     nTotalLG += nTotalSMLG[iSM]; 
381     nTotalHG += nTotalSMHG[iSM]; 
382     GetRawsData(kNsmodLG)->Fill(iSM, nTotalSMLG[iSM]); 
383     GetRawsData(kNsmodHG)->Fill(iSM, nTotalSMHG[iSM]); 
384   }
385   GetRawsData(kNtotLG)->Fill(nTotalLG);
386   GetRawsData(kNtotHG)->Fill(nTotalHG);
387
388   // just in case the next rawreader consumer forgets to reset; let's do it here again..
389   rawReader->Reset() ;
390
391   return;
392 }
393
394 //____________________________________________________________________________
395 void AliEMCALQADataMakerRec::MakeDigits(TClonesArray * digits)
396 {
397   // makes data from Digits
398  
399   // Check id histograms already created for this Event Specie
400   if ( ! GetDigitsData(0) )
401     InitDigits() ;
402
403   GetDigitsData(1)->Fill(digits->GetEntriesFast()) ; 
404   TIter next(digits) ; 
405   AliEMCALDigit * digit ; 
406   while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) {
407     GetDigitsData(0)->Fill( digit->GetAmp()) ;
408   }  
409   
410 }
411
412 //____________________________________________________________________________
413 void AliEMCALQADataMakerRec::MakeDigits(TTree * digitTree)
414 {
415   // makes data from Digit Tree
416   TClonesArray * digits = new TClonesArray("AliEMCALDigit", 1000) ; 
417   
418   TBranch * branch = digitTree->GetBranch("EMCAL") ;
419   if ( ! branch ) {
420     AliWarning("EMCAL branch in Digit Tree not found") ; 
421   } else {
422     branch->SetAddress(&digits) ;
423     branch->GetEntry(0) ; 
424     MakeDigits(digits) ; 
425   }
426   
427 }
428
429 //____________________________________________________________________________
430 void AliEMCALQADataMakerRec::MakeRecPoints(TTree * clustersTree)
431 {
432   // makes data from RecPoints
433   TBranch *emcbranch = clustersTree->GetBranch("EMCALECARP");
434   if (!emcbranch) { 
435     AliError("can't get the branch with the EMCAL clusters !");
436     return;
437   }
438
439   // Check id histograms already created for this Event Specie
440   if ( ! GetRecPointsData(kRecPM) )
441     InitRecPoints() ;
442   
443   TObjArray * emcrecpoints = new TObjArray(100) ;
444   emcbranch->SetAddress(&emcrecpoints);
445   emcbranch->GetEntry(0);
446   
447   GetRecPointsData(kRecPM)->Fill(emcrecpoints->GetEntriesFast()) ; 
448   TIter next(emcrecpoints) ; 
449   AliEMCALRecPoint * rp ; 
450   while ( (rp = dynamic_cast<AliEMCALRecPoint *>(next())) ) {
451     GetRecPointsData(kRecPE)->Fill( rp->GetEnergy()) ;
452     GetRecPointsData(kRecPDigM)->Fill(rp->GetMultiplicity());
453   }
454   emcrecpoints->Delete();
455   delete emcrecpoints;
456   
457 }
458
459 //____________________________________________________________________________ 
460 void AliEMCALQADataMakerRec::StartOfDetectorCycle()
461 {
462   //Detector specific actions at start of cycle
463   
464 }
465