use the new AliCaloRawStreamV3 decoder (from Yuri), based on the new AliAltroRawStrea...
[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 "AliCaloRawStreamV3.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   AliCaloRawStreamV3 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 maxTime = 0;
285   int startBin = 0;
286
287   // calc. quantities
288   double meanPed = 0, squaredMean = 0, rmsPed = 0;
289
290   // start loop over input stream  
291   while (in.NextDDL()) {
292     while (in.NextChannel()) {
293
294       // counters
295       int max = sampleMin, min = sampleMax; // min and max sample values
296       
297       // for the pedestal calculation
298       int sampleSum = 0; // sum of samples
299       int squaredSampleSum = 0; // sum of samples squared
300       int nSum = 0; // number of samples in sum
301       
302       while (in.NextBunch()) {
303         const UShort_t *sig = in.GetSignals();
304         startBin = in.GetStartTimeBin();
305         for (i = 0; i < in.GetBunchLength(); i++) {
306           sample = sig[i];
307           time = startBin--;
308
309           // check if it's a min or max value
310           if (sample < min) min = sample;
311           if (sample > max) {
312             max = sample;
313             maxTime = time;
314           }
315
316           // should we add it for the pedestal calculation?
317           if ( (firstPedestalSample<=time && time<=lastPedestalSample) || // sample time in range
318                !selectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all 
319             sampleSum += sample;
320             squaredSampleSum += sample*sample;
321             nSum++;
322           }
323           
324         } // loop over samples in bunch
325       } // loop over bunches
326
327       // calculate pedesstal estimate: mean of possibly selected samples
328       if (nSum > 0) {
329         meanPed = sampleSum / (1.0 * nSum);
330         squaredMean = squaredSampleSum / (1.0 * nSum);
331         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
332         rmsPed = sqrt(squaredMean - meanPed*meanPed); 
333       }
334       else {
335         meanPed = 0;
336         squaredMean = 0;
337         rmsPed  = 0;
338       }
339
340       //If we're here then we're done with this tower
341       gain = -1; // init to not valid value
342       if ( in.IsLowGain() ) {
343         gain = 0;
344       }
345       else if ( in.IsHighGain() ) {
346         gain = 1;
347       }
348       
349       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
350       iSM = in.GetModule(); //The modules are numbered starting from 0
351
352       if (iSM>=0 && iSM<fSuperModules) { // valid module reading, can go on with filling
353
354         int towerId = iSM*nTowersPerSM + in.GetColumn()*nRows + in.GetRow();
355
356         if (gain == 0) { 
357           //fill the low gain histograms, and counters
358           nTotalSMLG[iSM]++; // one more channel found
359           GetRawsData(kSigLG)->Fill(towerId, max - min);
360           GetRawsData(kTimeLG)->Fill(towerId, maxTime);
361           if (nSum>0) { // only fill pedestal info in case it could be calculated
362             GetRawsData(kPedLG)->Fill(towerId, meanPed);
363             GetRawsData(kPedRMSLG)->Fill(towerId, rmsPed);
364           }
365         } // gain==0
366         else if (gain == 1) {           
367           //fill the high gain ones
368           nTotalSMHG[iSM]++; // one more channel found
369           GetRawsData(kSigHG)->Fill(towerId, max - min);
370           GetRawsData(kTimeHG)->Fill(towerId, maxTime);
371           if (nSum>0) { // only fill pedestal info in case it could be calculated
372             GetRawsData(kPedHG)->Fill(towerId, meanPed);
373             GetRawsData(kPedRMSHG)->Fill(towerId, rmsPed);
374           }
375         }
376       } // SM index OK
377
378     }// end while over channel 
379    
380   }//end while over DDL's, of input stream 
381
382   // let's also fill the SM and event counter histograms
383   int nTotalHG = 0;
384   int nTotalLG = 0;
385   for (iSM=0; iSM<fSuperModules; iSM++) {  
386     nTotalLG += nTotalSMLG[iSM]; 
387     nTotalHG += nTotalSMHG[iSM]; 
388     GetRawsData(kNsmodLG)->Fill(iSM, nTotalSMLG[iSM]); 
389     GetRawsData(kNsmodHG)->Fill(iSM, nTotalSMHG[iSM]); 
390   }
391   GetRawsData(kNtotLG)->Fill(nTotalLG);
392   GetRawsData(kNtotHG)->Fill(nTotalHG);
393
394   // just in case the next rawreader consumer forgets to reset; let's do it here again..
395   rawReader->Reset() ;
396
397   return;
398 }
399
400 //____________________________________________________________________________
401 void AliEMCALQADataMakerRec::MakeDigits(TClonesArray * digits)
402 {
403   // makes data from Digits
404  
405   // Check id histograms already created for this Event Specie
406   if ( ! GetDigitsData(0) )
407     InitDigits() ;
408
409   GetDigitsData(1)->Fill(digits->GetEntriesFast()) ; 
410   TIter next(digits) ; 
411   AliEMCALDigit * digit ; 
412   while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) {
413     GetDigitsData(0)->Fill( digit->GetAmp()) ;
414   }  
415   
416 }
417
418 //____________________________________________________________________________
419 void AliEMCALQADataMakerRec::MakeDigits(TTree * digitTree)
420 {
421   // makes data from Digit Tree
422   TClonesArray * digits = new TClonesArray("AliEMCALDigit", 1000) ; 
423   
424   TBranch * branch = digitTree->GetBranch("EMCAL") ;
425   if ( ! branch ) {
426     AliWarning("EMCAL branch in Digit Tree not found") ; 
427   } else {
428     branch->SetAddress(&digits) ;
429     branch->GetEntry(0) ; 
430     MakeDigits(digits) ; 
431   }
432   
433 }
434
435 //____________________________________________________________________________
436 void AliEMCALQADataMakerRec::MakeRecPoints(TTree * clustersTree)
437 {
438   // makes data from RecPoints
439   TBranch *emcbranch = clustersTree->GetBranch("EMCALECARP");
440   if (!emcbranch) { 
441     AliError("can't get the branch with the EMCAL clusters !");
442     return;
443   }
444
445   // Check id histograms already created for this Event Specie
446   if ( ! GetRecPointsData(kRecPM) )
447     InitRecPoints() ;
448   
449   TObjArray * emcrecpoints = new TObjArray(100) ;
450   emcbranch->SetAddress(&emcrecpoints);
451   emcbranch->GetEntry(0);
452   
453   GetRecPointsData(kRecPM)->Fill(emcrecpoints->GetEntriesFast()) ; 
454   TIter next(emcrecpoints) ; 
455   AliEMCALRecPoint * rp ; 
456   while ( (rp = dynamic_cast<AliEMCALRecPoint *>(next())) ) {
457     GetRecPointsData(kRecPE)->Fill( rp->GetEnergy()) ;
458     GetRecPointsData(kRecPDigM)->Fill(rp->GetMultiplicity());
459   }
460   emcrecpoints->Delete();
461   delete emcrecpoints;
462   
463 }
464
465 //____________________________________________________________________________ 
466 void AliEMCALQADataMakerRec::StartOfDetectorCycle()
467 {
468   //Detector specific actions at start of cycle
469   
470 }
471