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