]> git.uio.no Git - u/mrichter/AliRoot.git/blob - EMCAL/AliEMCALQADataMakerRec.cxx
Warning corrected about fSuperModules not initialized in member initialization list
[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",    200, 0., 20.) ; 
94   h1->Sumw2() ;
95   Add2ESDsList(h1, kESDCaloClusE, !expert, image)  ;                                                     
96
97   TH1I * h2 = new TH1I("hESDCaloClusterM", "ESDs CaloCluster multiplicity in EMCAL", 100, 0,  100) ; 
98   h2->Sumw2() ;
99   Add2ESDsList(h2, kESDCaloClusM, !expert, image)  ;
100
101   TH1F * h3 = new TH1F("hESDCaloCellA",  "ESDs CaloCell amplitude in EMCAL",    500, 0., 250.) ; 
102   h3->Sumw2() ;
103   Add2ESDsList(h3, kESDCaloCellA, !expert, image)  ;  
104  
105   TH1I * h4 = new TH1I("hESDCaloCellM", "ESDs CaloCell multiplicity in EMCAL", 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",    500, 0, 500) ; 
119   h0->Sumw2() ;
120   Add2DigitsList(h0, 0, !expert, image) ;
121   TH1I * h1 = new TH1I("hEmcalDigitsMul", "Digits multiplicity distribution in EMCAL", 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",200, 0.,20.); //GeV
134   h0->Sumw2();
135   Add2RecPointsList(h0,kRecPE, !expert, image);
136
137   TH1I* h1 = new TH1I("hEMCALRpM","EMCAL RecPoint multiplicities",100,0,100);
138   h1->Sumw2();
139   Add2RecPointsList(h1,kRecPM, !expert, image);
140
141   TH1I* h2 = new TH1I("hEMCALRpDigM","EMCAL RecPoint Digit Multiplicities",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",
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",  
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", 
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", 
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", 
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",
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", 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", 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", 
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",
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", 
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",
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   Int_t nTot = 0 ; 
214   for ( Int_t index = 0; index < esd->GetNumberOfCaloClusters() ; index++ ) {
215     AliESDCaloCluster * clu = esd->GetCaloCluster(index) ;
216     if( clu->IsEMCAL() ) {
217       GetESDsData(kESDCaloClusE)->Fill(clu->E()) ;
218       nTot++ ;
219     } 
220   }
221   GetESDsData(kESDCaloClusM)->Fill(nTot) ;
222
223   //fill calo cells
224   AliESDCaloCells* cells = esd->GetEMCALCells();
225   GetESDsData(kESDCaloCellM)->Fill(cells->GetNumberOfCells()) ;
226
227   for ( Int_t index = 0; index < cells->GetNumberOfCells() ; index++ ) {
228     GetESDsData(kESDCaloCellA)->Fill(cells->GetAmplitude(index)) ;
229   }
230
231 }
232
233 //____________________________________________________________________________
234 void AliEMCALQADataMakerRec::MakeRaws(AliRawReader* rawReader)
235 {
236   //Fill prepared histograms with Raw digit properties
237
238   //Raw histogram filling not yet implemented
239   //
240   //Need to figure out how to get the info we want without having to
241   //actually run Raw2Digits twice.
242   //I suspect what we actually want is a raw digits method, not a true
243   //emcal raw data method, but this doesn't seem to be allowed in
244   //AliQADataMakerRec.h
245
246   // For now, to avoid redoing the expensive signal fits we just
247   // look at max vs min of the signal spextra, a la online usage in
248   // AliCaloCalibPedestal
249
250   rawReader->Reset() ;
251   AliCaloRawStream in(rawReader,"EMCAL"); 
252
253   // setup
254   int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
255   int nRows = 24; // number of rows per SuperModule
256   int sampleMin = 0; 
257   int sampleMax = 0x3ff; // 1023 = 10-bit range
258
259   // for the pedestal calculation
260   Bool_t selectPedestalSamples = kTRUE;
261   int firstPedestalSample = 0;
262   int lastPedestalSample = 15;
263
264   // SM counters; decl. should be safe, assuming we don't get more than 12 SuperModules..
265   int nTotalSMLG[12] = {0};
266   int nTotalSMHG[12] = {0};
267
268   // indices for the reading
269   int iSM = 0;
270   int sample = 0;
271   int gain = 0;
272   int time = 0;
273   // counters, on sample level
274   int i = 0; // the sample number in current event.
275   int max = sampleMin, min = sampleMax;//Use these for picking the pedestal
276   int maxTime = 0;
277
278   // for the pedestal calculation
279   int sampleSum = 0; // sum of samples
280   int squaredSampleSum = 0; // sum of samples squared
281   int nSum = 0; // number of samples in sum
282   // calc. quantities
283   double meanPed = 0, squaredMean = 0, rmsPed = 0;
284   
285   while (in.Next()) { // loop over input stream
286     sample = in.GetSignal(); //Get the adc signal
287     time = in.GetTime();
288     if (sample < min) { min = sample; }
289     if (sample > max) { 
290       max = sample;
291       maxTime = time;
292     }
293     i++;
294
295     // should we add it for the pedestal calculation?
296     if ( (firstPedestalSample<=time && time<=lastPedestalSample) || // sample time in range
297          !selectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all 
298       sampleSum += sample;
299       squaredSampleSum += sample*sample;
300       nSum++;
301     }
302
303     if ( i >= in.GetTimeLength()) {
304       // calculate pedesstal estimate: mean of possibly selected samples
305       if (nSum > 0) {
306         meanPed = sampleSum / (1.0 * nSum);
307         squaredMean = squaredSampleSum / (1.0 * nSum);
308         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
309         rmsPed = sqrt(squaredMean - meanPed*meanPed); 
310       }
311       else {
312         meanPed = 0;
313         squaredMean = 0;
314         rmsPed  = 0;
315       }
316
317       //If we're here then we're done with this tower
318       gain = -1; // init to not valid value
319       if ( in.IsLowGain() ) {
320         gain = 0;
321       }
322       else if ( in.IsHighGain() ) {
323         gain = 1;
324       }
325       
326       iSM = in.GetModule(); //The modules are numbered starting from 0
327
328       if (iSM>=0 && iSM<fSuperModules) { // valid module reading, can go on with filling
329
330         int towerId = iSM*nTowersPerSM + in.GetColumn()*nRows + in.GetRow();
331
332         if (gain == 0) { 
333           //fill the low gain histograms, and counters
334           nTotalSMLG[iSM]++; // one more channel found
335           GetRawsData(kSigLG)->Fill(towerId, max - min);
336           GetRawsData(kTimeLG)->Fill(towerId, maxTime);
337           if (nSum>0) { // only fill pedestal info in case it could be calculated
338             GetRawsData(kPedLG)->Fill(towerId, meanPed);
339             GetRawsData(kPedRMSLG)->Fill(towerId, rmsPed);
340           }
341         } // gain==0
342         else if (gain == 1) {           
343           //fill the high gain ones
344           nTotalSMHG[iSM]++; // one more channel found
345           GetRawsData(kSigHG)->Fill(towerId, max - min);
346           GetRawsData(kTimeHG)->Fill(towerId, maxTime);
347           if (nSum>0) { // only fill pedestal info in case it could be calculated
348             GetRawsData(kPedHG)->Fill(towerId, meanPed);
349             GetRawsData(kPedRMSHG)->Fill(towerId, rmsPed);
350           }
351         }
352       } // SM index OK
353
354       // reset counters
355       max = sampleMin; min = sampleMax;
356       maxTime = 0;
357       i = 0;
358       // also pedestal calc counters
359       sampleSum = 0; // sum of samples
360       squaredSampleSum = 0; // sum of samples squared
361       nSum = 0; // number of samples in sum
362     
363     }//End if, of channel
364    
365   }//end while, of stream
366
367   // let's also fill the SM and event counter histograms
368   int nTotalHG = 0;
369   int nTotalLG = 0;
370   for (iSM=0; iSM<fSuperModules; iSM++) {  
371     nTotalLG += nTotalSMLG[iSM]; 
372     nTotalHG += nTotalSMHG[iSM]; 
373     GetRawsData(kNsmodLG)->Fill(iSM, nTotalSMLG[iSM]); 
374     GetRawsData(kNsmodHG)->Fill(iSM, nTotalSMHG[iSM]); 
375   }
376   GetRawsData(kNtotLG)->Fill(nTotalLG);
377   GetRawsData(kNtotHG)->Fill(nTotalHG);
378
379   // just in case the next rawreader consumer forgets to reset; let's do it here again..
380   rawReader->Reset() ;
381
382   return;
383 }
384
385 //____________________________________________________________________________
386 void AliEMCALQADataMakerRec::MakeDigits(TClonesArray * digits)
387 {
388   // makes data from Digits
389   
390   GetDigitsData(1)->Fill(digits->GetEntriesFast()) ; 
391   TIter next(digits) ; 
392   AliEMCALDigit * digit ; 
393   while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) {
394     GetDigitsData(0)->Fill( digit->GetAmp()) ;
395   }  
396   
397 }
398
399 //____________________________________________________________________________
400 void AliEMCALQADataMakerRec::MakeDigits(TTree * digitTree)
401 {
402   // makes data from Digit Tree
403   TClonesArray * digits = new TClonesArray("AliEMCALDigit", 1000) ; 
404   
405   TBranch * branch = digitTree->GetBranch("EMCAL") ;
406   if ( ! branch ) {
407     AliWarning("EMCAL branch in Digit Tree not found") ; 
408   } else {
409     branch->SetAddress(&digits) ;
410     branch->GetEntry(0) ; 
411     MakeDigits(digits) ; 
412   }
413   
414 }
415
416 //____________________________________________________________________________
417 void AliEMCALQADataMakerRec::MakeRecPoints(TTree * clustersTree)
418 {
419   // makes data from RecPoints
420   TBranch *emcbranch = clustersTree->GetBranch("EMCALECARP");
421   if (!emcbranch) { 
422     AliError("can't get the branch with the EMCAL clusters !");
423     return;
424   }
425   TObjArray * emcrecpoints = new TObjArray(100) ;
426   emcbranch->SetAddress(&emcrecpoints);
427   emcbranch->GetEntry(0);
428   
429   GetRecPointsData(kRecPM)->Fill(emcrecpoints->GetEntriesFast()) ; 
430   TIter next(emcrecpoints) ; 
431   AliEMCALRecPoint * rp ; 
432   while ( (rp = dynamic_cast<AliEMCALRecPoint *>(next())) ) {
433     GetRecPointsData(kRecPE)->Fill( rp->GetEnergy()) ;
434     GetRecPointsData(kRecPDigM)->Fill(rp->GetMultiplicity());
435   }
436   emcrecpoints->Delete();
437   delete emcrecpoints;
438   
439 }
440
441 //____________________________________________________________________________ 
442 void AliEMCALQADataMakerRec::StartOfDetectorCycle()
443 {
444   //Detector specific actions at start of cycle
445   
446 }
447