]>
Commit | Line | Data |
---|---|---|
94594e5d | 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> | |
3fe36ccb | 29 | #include <TProfile.h> |
94594e5d | 30 | |
31 | // --- Standard library --- | |
32 | ||
3fe36ccb | 33 | |
94594e5d | 34 | // --- AliRoot header files --- |
35 | #include "AliESDCaloCluster.h" | |
601c73e3 | 36 | #include "AliESDCaloCells.h" |
94594e5d | 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" | |
78328afd | 45 | #include "AliRawReader.h" |
3fe36ccb | 46 | #include "AliCaloRawStream.h" |
94594e5d | 47 | |
48 | ClassImp(AliEMCALQADataMakerRec) | |
49 | ||
50 | //____________________________________________________________________________ | |
51 | AliEMCALQADataMakerRec::AliEMCALQADataMakerRec() : | |
3fe36ccb | 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.. | |
94594e5d | 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()); | |
3fe36ccb | 65 | fSuperModules = qadm.GetSuperModules(); |
94594e5d | 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 | //____________________________________________________________________________ | |
4e25ac79 | 78 | void AliEMCALQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list) |
94594e5d | 79 | { |
80 | //Detector specific actions at end of cycle | |
81 | // do the QA checking | |
4e25ac79 | 82 | AliQAChecker::Instance()->Run(AliQAv1::kEMCAL, task, list) ; |
94594e5d | 83 | } |
84 | ||
85 | //____________________________________________________________________________ | |
86 | void AliEMCALQADataMakerRec::InitESDs() | |
87 | { | |
88 | //Create histograms to controll ESD | |
7d297381 | 89 | const Bool_t expert = kTRUE ; |
90 | const Bool_t image = kTRUE ; | |
91 | ||
601c73e3 | 92 | TH1F * h1 = new TH1F("hESDCaloClusterE", "ESDs CaloCluster energy in EMCAL", 200, 0., 20.) ; |
94594e5d | 93 | h1->Sumw2() ; |
7d297381 | 94 | Add2ESDsList(h1, kESDCaloClusE, !expert, image) ; |
601c73e3 | 95 | |
96 | TH1I * h2 = new TH1I("hESDCaloClusterM", "ESDs CaloCluster multiplicity in EMCAL", 100, 0, 100) ; | |
94594e5d | 97 | h2->Sumw2() ; |
7d297381 | 98 | Add2ESDsList(h2, kESDCaloClusM, !expert, image) ; |
601c73e3 | 99 | |
100 | TH1F * h3 = new TH1F("hESDCaloCellA", "ESDs CaloCell amplitude in EMCAL", 500, 0., 250.) ; | |
94594e5d | 101 | h3->Sumw2() ; |
7d297381 | 102 | Add2ESDsList(h3, kESDCaloCellA, !expert, image) ; |
94594e5d | 103 | |
601c73e3 | 104 | TH1I * h4 = new TH1I("hESDCaloCellM", "ESDs CaloCell multiplicity in EMCAL", 200, 0, 1000) ; |
94594e5d | 105 | h4->Sumw2() ; |
7d297381 | 106 | Add2ESDsList(h4, kESDCaloCellM, !expert, image) ; |
94594e5d | 107 | |
108 | } | |
109 | ||
110 | //____________________________________________________________________________ | |
111 | void AliEMCALQADataMakerRec::InitRecPoints() | |
112 | { | |
113 | // create Reconstructed Points histograms in RecPoints subdir | |
7d297381 | 114 | const Bool_t expert = kTRUE ; |
115 | const Bool_t image = kTRUE ; | |
116 | ||
601c73e3 | 117 | TH1F* h0 = new TH1F("hEMCALRpE","EMCAL RecPoint energies",200, 0.,20.); //GeV |
118 | h0->Sumw2(); | |
7d297381 | 119 | Add2RecPointsList(h0,kRecPE, !expert, image); |
94594e5d | 120 | |
601c73e3 | 121 | TH1I* h1 = new TH1I("hEMCALRpM","EMCAL RecPoint multiplicities",100,0,100); |
122 | h1->Sumw2(); | |
7d297381 | 123 | Add2RecPointsList(h1,kRecPM, !expert, image); |
94594e5d | 124 | |
601c73e3 | 125 | TH1I* h2 = new TH1I("hEMCALRpDigM","EMCAL RecPoint Digit Multiplicities",20,0,20); |
126 | h2->Sumw2(); | |
7d297381 | 127 | Add2RecPointsList(h2,kRecPDigM, !expert, image); |
94594e5d | 128 | |
129 | } | |
130 | ||
131 | //____________________________________________________________________________ | |
132 | void AliEMCALQADataMakerRec::InitRaws() | |
133 | { | |
134 | // create Raws histograms in Raws subdir | |
7d297381 | 135 | const Bool_t expert = kTRUE ; |
136 | const Bool_t saveCorr = kTRUE ; | |
137 | const Bool_t image = kTRUE ; | |
3fe36ccb | 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) ; | |
7d297381 | 145 | Add2RawsList(h0, kNsmodLG, !expert, image, !saveCorr) ; |
3fe36ccb | 146 | TProfile * h1 = new TProfile("hHighEmcalSupermodules", "High Gain EMC: # of towers vs SuperMod", |
147 | fSuperModules, -0.5, fSuperModules-0.5) ; | |
7d297381 | 148 | Add2RawsList(h1, kNsmodHG, !expert, image, !saveCorr) ; |
94594e5d | 149 | |
3fe36ccb | 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) ; | |
94594e5d | 157 | |
3fe36ccb | 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) ; | |
601c73e3 | 165 | |
3fe36ccb | 166 | // total counter: channels per event |
167 | TH1I * h6 = new TH1I("hLowNtot", "Low Gain EMC: Total Number of found towers", 200, 0, nTot) ; | |
601c73e3 | 168 | h6->Sumw2() ; |
7d297381 | 169 | Add2RawsList(h6, kNtotLG, !expert, image, !saveCorr) ; |
3fe36ccb | 170 | TH1I * h7 = new TH1I("hHighNtot", "High Gain EMC: Total Number of found towers", 200,0, nTot) ; |
601c73e3 | 171 | h7->Sumw2() ; |
7d297381 | 172 | Add2RawsList(h7, kNtotHG, !expert, image, !saveCorr) ; |
601c73e3 | 173 | |
3fe36ccb | 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) ; | |
94594e5d | 189 | |
190 | } | |
191 | ||
192 | //____________________________________________________________________________ | |
193 | void AliEMCALQADataMakerRec::MakeESDs(AliESDEvent * esd) | |
194 | { | |
195 | // make QA data from ESDs | |
196 | ||
197 | Int_t nTot = 0 ; | |
94594e5d | 198 | for ( Int_t index = 0; index < esd->GetNumberOfCaloClusters() ; index++ ) { |
199 | AliESDCaloCluster * clu = esd->GetCaloCluster(index) ; | |
200 | if( clu->IsEMCAL() ) { | |
601c73e3 | 201 | GetESDsData(kESDCaloClusE)->Fill(clu->E()) ; |
94594e5d | 202 | nTot++ ; |
203 | } | |
204 | } | |
601c73e3 | 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 | ||
94594e5d | 215 | } |
216 | ||
217 | //____________________________________________________________________________ | |
78328afd | 218 | void AliEMCALQADataMakerRec::MakeRaws(AliRawReader* rawReader) |
94594e5d | 219 | { |
220 | //Fill prepared histograms with Raw digit properties | |
221 | ||
222 | //Raw histogram filling not yet implemented | |
601c73e3 | 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 | |
94594e5d | 229 | |
3fe36ccb | 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; | |
94594e5d | 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 | ||
601c73e3 | 382 | GetRecPointsData(kRecPM)->Fill(emcrecpoints->GetEntriesFast()) ; |
94594e5d | 383 | TIter next(emcrecpoints) ; |
384 | AliEMCALRecPoint * rp ; | |
94594e5d | 385 | while ( (rp = dynamic_cast<AliEMCALRecPoint *>(next())) ) { |
601c73e3 | 386 | GetRecPointsData(kRecPE)->Fill( rp->GetEnergy()) ; |
387 | GetRecPointsData(kRecPDigM)->Fill(rp->GetMultiplicity()); | |
94594e5d | 388 | } |
94594e5d | 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 | } | |
3fe36ccb | 400 |