1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //________________________________________________________________________
19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
20 // It can be created and used a la (ctor):
22 //Create the object for making the histograms
23 fPedestals = new AliCaloCalibPedestal( fDetType );
24 // AliCaloCalibPedestal knows how many modules we have for PHOS or EMCAL
25 fNumModules = fPedestals->GetModules();
28 // fPedestals->ProcessEvent(fCaloRawStream);
29 // asked to draw histograms:
30 // fPedestals->GetDeadMap(i)->Draw("col");
32 // fPedestals->GetPeakProfileHighGainRatio((i < fNumModules) ? i : fVisibleModule)->Draw("colz");
34 // The pseudo-code examples above were from the first implementation in MOOD (summer 2007).
35 //________________________________________________________________________
37 //#include "TCanvas.h"
47 #include "AliRawReader.h"
48 #include "AliCaloRawStreamV3.h"
51 #include "AliCaloCalibPedestal.h"
53 ClassImp(AliCaloCalibPedestal)
57 // ctor; initialize everything in order to avoid compiler warnings
58 AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
62 fPedestalLEDRefLowGain(),
63 fPedestalLEDRefHighGain(),
64 fPeakMinusPedLowGain(),
65 fPeakMinusPedHighGain(),
66 fPeakMinusPedHighGainHisto(),
67 fPedestalLowGainDiff(),
68 fPedestalHighGainDiff(),
69 fPedestalLEDRefLowGainDiff(),
70 fPedestalLEDRefHighGainDiff(),
71 fPeakMinusPedLowGainDiff(),
72 fPeakMinusPedHighGainDiff(),
73 fPedestalLowGainRatio(),
74 fPedestalHighGainRatio(),
75 fPedestalLEDRefLowGainRatio(),
76 fPedestalLEDRefHighGainRatio(),
77 fPeakMinusPedLowGainRatio(),
78 fPeakMinusPedHighGainRatio(),
84 fResurrectedTowers(0),
97 fSelectPedestalSamples(kTRUE),
98 fFirstPedestalSample(0),
99 fLastPedestalSample(15),
101 fWarningThreshold(50),
102 fWarningFraction(0.002),
105 //Default constructor. First we set the detector-type related constants.
106 if (detectorType == kPhos) {
107 fColumns = fgkPhosCols;
109 fLEDRefs = fgkPhosLEDRefs;
110 fModules = fgkPhosModules;
111 fCaloString = "PHOS";
117 //We'll just trust the enum to keep everything in line, so that if detectorType
118 //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
119 //case, if someone intentionally gives another number
120 fColumns = AliEMCALGeoParams::fgkEMCALCols;
121 fRows = AliEMCALGeoParams::fgkEMCALRows;
122 fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
123 fModules = AliEMCALGeoParams::fgkEMCALModules;
124 fCaloString = "EMCAL";
129 fDetType = detectorType;
131 //Then, loop for the requested number of modules
133 for (int i = 0; i < fModules; i++) {
134 //Pedestals, low gain
135 name = "hPedlowgain";
137 title = "Pedestals, low gain, module ";
139 fPedestalLowGain.Add(new TProfile2D(name, title,
140 fColumns, 0.0, fColumns,
141 fRows, fRowMin, fRowMax,"s"));
143 //Pedestals, high gain
144 name = "hPedhighgain";
146 title = "Pedestals, high gain, module ";
148 fPedestalHighGain.Add(new TProfile2D(name, title,
149 fColumns, 0.0, fColumns,
150 fRows, fRowMin, fRowMax,"s"));
152 //LED Ref/Mon pedestals, low gain
153 name = "hPedestalLEDReflowgain";
155 title = "Pedestal LEDRef, low gain, module ";
157 fPedestalLEDRefLowGain.Add(new TProfile(name, title,
158 fLEDRefs, 0.0, fLEDRefs, "s"));
160 //LED Ref/Mon pedestals, high gain
161 name = "hPedestalLEDRefhighgain";
163 title = "Pedestal LEDRef, high gain, module ";
165 fPedestalLEDRefHighGain.Add(new TProfile(name, title,
166 fLEDRefs, 0.0, fLEDRefs, "s"));
168 //Peak-Pedestals, low gain
169 name = "hPeakMinusPedlowgain";
171 title = "Peak-Pedestal, low gain, module ";
173 fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
174 fColumns, 0.0, fColumns,
175 fRows, fRowMin, fRowMax,"s"));
177 //Peak-Pedestals, high gain
178 name = "hPeakMinusPedhighgain";
180 title = "Peak-Pedestal, high gain, module ";
182 fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
183 fColumns, 0.0, fColumns,
184 fRows, fRowMin, fRowMax,"s"));
186 //Peak-Pedestals, high gain - TH2F histo
187 name = "hPeakMinusPedhighgainHisto";
189 title = "Peak-Pedestal, high gain, module ";
191 fPeakMinusPedHighGainHisto.Add(new TH2F(name, title,
192 fColumns*fRows, 0.0, fColumns*fRows,
197 title = "Dead map, module ";
199 fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
200 fRows, fRowMin, fRowMax));
202 }//end for nModules create the histograms
204 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
205 fPedestalLowGain.Compress();
206 fPedestalHighGain.Compress();
207 fPedestalLEDRefLowGain.Compress();
208 fPedestalLEDRefHighGain.Compress();
209 fPeakMinusPedLowGain.Compress();
210 fPeakMinusPedHighGain.Compress();
211 fPeakMinusPedHighGainHisto.Compress();
217 //_____________________________________________________________________
218 AliCaloCalibPedestal::~AliCaloCalibPedestal()
220 if (fReference) delete fReference;//Delete the reference object, if it has been loaded
221 //TObjArray will delete the histos/profiles when it is deleted.
225 //_____________________________________________________________________
226 AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
230 fPedestalLEDRefLowGain(),
231 fPedestalLEDRefHighGain(),
232 fPeakMinusPedLowGain(),
233 fPeakMinusPedHighGain(),
234 fPeakMinusPedHighGainHisto(),
235 fPedestalLowGainDiff(),
236 fPedestalHighGainDiff(),
237 fPedestalLEDRefLowGainDiff(),
238 fPedestalLEDRefHighGainDiff(),
239 fPeakMinusPedLowGainDiff(),
240 fPeakMinusPedHighGainDiff(),
241 fPedestalLowGainRatio(),
242 fPedestalHighGainRatio(),
243 fPedestalLEDRefLowGainRatio(),
244 fPedestalLEDRefHighGainRatio(),
245 fPeakMinusPedLowGainRatio(),
246 fPeakMinusPedHighGainRatio(),
248 fNEvents(ped.GetNEvents()),
249 fNChanFills(ped.GetNChanFills()),
250 fDeadTowers(ped.GetDeadTowerCount()),
251 fNewDeadTowers(ped.GetDeadTowerNew()),
252 fResurrectedTowers(ped.GetDeadTowerResurrected()),
253 fReference( 0 ), //! note that we do not try to copy the reference info here
254 fDetType(ped.GetDetectorType()),
255 fColumns(ped.GetColumns()),
256 fRows(ped.GetRows()),
257 fLEDRefs(ped.GetLEDRefs()),
258 fModules(ped.GetModules()),
259 fRowMin(ped.GetRowMin()),
260 fRowMax(ped.GetRowMax()),
261 fRowMultiplier(ped.GetRowMultiplier()),
262 fCaloString(ped.GetCaloString()),
263 fMapping(NULL), //! note that we are not copying the map info
264 fRunNumber(ped.GetRunNumber()),
265 fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
266 fFirstPedestalSample(ped.GetFirstPedestalSample()),
267 fLastPedestalSample(ped.GetLastPedestalSample()),
268 fDeadThreshold(ped.GetDeadThreshold()),
269 fWarningThreshold(ped.GetWarningThreshold()),
270 fWarningFraction(ped.GetWarningFraction()),
271 fHotSigma(ped.GetHotSigma())
273 // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
274 //DS: this has not really been tested yet..
275 for (int i = 0; i < fModules; i++) {
276 fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
277 fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
278 fPedestalLEDRefLowGain.Add( ped.GetPedLEDRefProfileLowGain(i) );
279 fPedestalLEDRefHighGain.Add( ped.GetPedLEDRefProfileHighGain(i) );
280 fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
281 fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
282 fPeakMinusPedHighGainHisto.Add( ped.GetPeakHighGainHisto(i) );
284 fDeadMap.Add( ped.GetDeadMap(i) );
287 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
288 fPedestalLowGain.Compress();
289 fPedestalHighGain.Compress();
290 fPedestalLEDRefLowGain.Compress();
291 fPedestalLEDRefHighGain.Compress();
292 fPeakMinusPedLowGain.Compress();
293 fPeakMinusPedHighGain.Compress();
294 fPeakMinusPedHighGainHisto.Compress();
299 // assignment operator; use copy ctor to make life easy..
300 //_____________________________________________________________________
301 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
303 // assignment operator; use copy ctor
304 if (&source == this) return *this;
306 new (this) AliCaloCalibPedestal(source);
310 //_____________________________________________________________________
311 void AliCaloCalibPedestal::Reset()
313 // Reset all arrays/histograms
314 for (int i = 0; i < fModules; i++) {
315 GetPedProfileLowGain(i)->Reset();
316 GetPedProfileHighGain(i)->Reset();
317 GetPedLEDRefProfileLowGain(i)->Reset();
318 GetPedLEDRefProfileHighGain(i)->Reset();
319 GetPeakProfileLowGain(i)->Reset();
320 GetPeakProfileHighGain(i)->Reset();
321 GetPeakHighGainHisto(i)->Reset();
322 GetDeadMap(i)->Reset();
324 if (!fPedestalLowGainDiff.IsEmpty()) {
325 //This means that the comparison profiles have been created.
327 GetPedProfileLowGainDiff(i)->Reset();
328 GetPedProfileHighGainDiff(i)->Reset();
329 GetPedLEDRefProfileLowGainDiff(i)->Reset();
330 GetPedLEDRefProfileHighGainDiff(i)->Reset();
331 GetPeakProfileLowGainDiff(i)->Reset();
332 GetPeakProfileHighGainDiff(i)->Reset();
334 GetPedProfileLowGainRatio(i)->Reset();
335 GetPedProfileHighGainRatio(i)->Reset();
336 GetPedLEDRefProfileLowGainRatio(i)->Reset();
337 GetPedLEDRefProfileHighGainRatio(i)->Reset();
338 GetPeakProfileLowGainRatio(i)->Reset();
339 GetPeakProfileHighGainRatio(i)->Reset();
346 fResurrectedTowers = 0;
348 //To think about: should fReference be deleted too?... let's not do it this time, at least...
351 // Parameter/cut handling
352 //_____________________________________________________________________
353 void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile)
355 // Note: this method is a bit more complicated than it really has to be
356 // - allowing for multiple entries per line, arbitrary order of the
357 // different variables etc. But I wanted to try and do this in as
358 // correct a C++ way as I could (as an exercise).
360 static const string delimitor("::");
362 // open, check input file
363 ifstream in( parameterFile );
365 printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
372 while ((in.rdstate() & ios::failbit) == 0 ) {
374 // Read into the raw char array and then construct a string
375 // to do the searching
376 in.getline(readline, 1024);
377 istringstream s(readline);
379 while ( ( s.rdstate() & ios::failbit ) == 0 ) {
384 // check stream status
385 if( s.rdstate() & ios::failbit ) break;
387 // skip rest of line if comments found
388 if( keyValue.substr( 0, 2 ) == "//" ) break;
390 // look for "::" in keyValue pair
391 size_t position = keyValue.find( delimitor );
392 if( position == string::npos ) {
393 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
396 // split keyValue pair
397 string key( keyValue.substr( 0, position ) );
398 string value( keyValue.substr( position+delimitor.size(),
399 keyValue.size()-delimitor.size() ) );
401 // check value does not contain a new delimitor
402 if( value.find( delimitor ) != string::npos ) {
403 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
406 // debug: check key value pair
407 // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
409 // if the key matches with something we expect, we assign the new value
410 istringstream iss(value);
411 // the comparison strings defined at the beginning of this method
412 if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") || (key == "fDeadThreshold") || (key == "fWarningThreshold") || (key == "fWarningFraction") || (key == "fHotSigma") ) {
413 printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
415 if (key == "fFirstPedestalSample") {
416 iss >> fFirstPedestalSample;
418 else if (key == "fLastPedestalSample") {
419 iss >> fLastPedestalSample;
421 else if (key == "fDeadThreshold") {
422 iss >> fDeadThreshold;
424 else if (key == "fWarningThreshold") {
425 iss >> fWarningThreshold;
427 else if (key == "fWarningFraction") {
428 iss >> fWarningFraction;
430 else if (key == "fHotSigma") {
444 //_____________________________________________________________________
445 void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
447 //Write parameters in file.
449 static const string delimitor("::");
450 ofstream out( parameterFile );
451 out << "// " << parameterFile << endl;
452 out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
453 out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
454 out << "fDeadThreshold" << "::" << fDeadThreshold << endl;
455 out << "fWarningThreshold" << "::" << fWarningThreshold << endl;
456 out << "fWarningFraction" << "::" << fWarningFraction << endl;
457 out << "fHotSigma" << "::" << fHotSigma << endl;
463 //_____________________________________________________________________
464 Bool_t AliCaloCalibPedestal::AddInfo(const AliCaloCalibPedestal *ped)
466 // just do this for the basic histograms/profiles that get filled in ProcessEvent
467 // may not have data for all modules, but let's just Add everything..
468 for (int i = 0; i < fModules; i++) {
469 GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
470 GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
471 GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
472 GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
473 GetPeakHighGainHisto(i)->Add( ped->GetPeakHighGainHisto(i) );
477 // DeadMap; Diff profiles etc would need to be redone after this operation
479 return kTRUE;//We succesfully added info from the supplied object
482 //_____________________________________________________________________
483 Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
485 // if fMapping is NULL the rawstream will crate its own mapping
486 AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
487 if (fDetType == kEmCal) {
488 rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
490 return ProcessEvent(&rawStream);
493 //_____________________________________________________________________
494 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
496 // Method to process=analyze one event in the data stream
497 if (!in) return kFALSE; //Return right away if there's a null pointer
498 fNEvents++; // one more event
500 // indices for the reading
503 int i = 0; // sample counter
506 // start loop over input stream
507 while (in->NextDDL()) {
508 while (in->NextChannel()) {
511 int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
516 vector<int> pedSamples;
518 while (in->NextBunch()) {
519 const UShort_t *sig = in->GetSignals();
520 startBin = in->GetStartTimeBin();
521 nsamples += in->GetBunchLength();
522 for (i = 0; i < in->GetBunchLength(); i++) {
526 // check if it's a min or max value
527 if (sample < min) min = sample;
528 if (sample > max) max = sample;
530 // should we add it for the pedestal calculation?
531 if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
532 !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
533 pedSamples.push_back( sig[i] );
537 } // loop over samples in bunch
538 } // loop over bunches
540 if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
542 // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
543 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
544 if (arrayPos >= fModules) {
545 //TODO: return an error message, if appopriate (perhaps if debug>0?)
549 if (arrayPos < 0 || arrayPos >= fModules) {
550 printf("Oh no: arrayPos = %i.\n", arrayPos);
553 fNChanFills++; // one more channel found, and profile to be filled
554 //NOTE: coordinates are (column, row) for the profiles
555 if ( in->IsLowGain() ) {
556 //fill the low gain histograms
557 ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
558 if (nPed>0) { // only fill pedestal info in case it could be calculated
559 for ( i=0; i<nPed; i++) {
560 ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]);
564 else if ( in->IsHighGain() ) {
565 //fill the high gain ones
566 ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
567 if (nPed>0) { // only fill pedestal info in case it could be calculated
568 for ( i=0; i<nPed; i++) {
569 ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]);
572 // for warning checks
573 int idx = in->GetRow() + fRows * in->GetColumn();
574 ((TH2F*)fPeakMinusPedHighGainHisto[arrayPos])->Fill(idx, max - min);
576 else if ( in->IsLEDMonData() ) {
577 // for LED Mon data, the mapping class holds the gain info in the Row variable
578 // and the Strip number in the Column..
579 int gain = in->GetRow();
580 int stripId = in->GetColumn();
581 if (nPed>0 && stripId<fLEDRefs) {
583 for ( i=0; i<nPed; i++) {
584 ((TProfile*)fPedestalLEDRefLowGain[arrayPos])->Fill(stripId, pedSamples[i]);
588 for ( i=0; i<nPed; i++) {
589 ((TProfile*)fPedestalLEDRefHighGain[arrayPos])->Fill(stripId, pedSamples[i]);
595 } // nsamples>0 check, some data found for this channel; not only trailer/header
596 }// end while over channel
597 }//end while over DDL's, of input stream
599 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
604 //_____________________________________________________________________
605 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
607 //Saves all the histograms (or profiles, to be accurate) to the designated file
609 TFile destFile(fileName, "recreate");
611 if (destFile.IsZombie()) {
617 for (int i = 0; i < fModules; i++) {
618 if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
619 fPeakMinusPedLowGain[i]->Write();
621 if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
622 fPeakMinusPedHighGain[i]->Write();
624 if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
625 fPedestalLowGain[i]->Write();
627 if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
628 fPedestalHighGain[i]->Write();
630 if( ((TProfile *)fPedestalLEDRefLowGain[i])->GetEntries() || saveEmptyHistos) {
631 fPedestalLEDRefLowGain[i]->Write();
633 if( ((TProfile *)fPedestalLEDRefHighGain[i])->GetEntries() || saveEmptyHistos) {
634 fPedestalLEDRefHighGain[i]->Write();
636 if( ((TH2F *)fPeakMinusPedHighGainHisto[i])->GetEntries() || saveEmptyHistos) {
637 fPeakMinusPedHighGainHisto[i]->Write();
647 //_____________________________________________________________________
648 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
651 //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
652 TH1::AddDirectory(kFALSE);
654 TFile *sourceFile = new TFile(fileName);
655 if (sourceFile->IsZombie()) {
656 return kFALSE;//We couldn't load the reference
659 if (fReference) delete fReference;//Delete the reference object, if it already exists
662 fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
664 if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
665 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
672 //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
673 //if we are called by someone who has set it to false...
674 TH1::AddDirectory(kTRUE);
676 return kTRUE;//We succesfully loaded the object
680 //_____________________________________________________________________
681 Bool_t AliCaloCalibPedestal::SetReference(AliCaloCalibPedestal *ref)
683 if (fReference) delete fReference;//Delete the reference object, if it already exists
688 if (!fReference || (fReference->GetDetectorType() != fDetType)) {
689 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
694 return kTRUE;//We succesfully loaded the object
697 //_____________________________________________________________________
698 void AliCaloCalibPedestal::ValidateComparisonProfiles()
700 //Make sure the comparison histos exist
701 if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
705 //Then, loop for the requested number of modules
707 for (int i = 0; i < fModules; i++) {
708 //Pedestals, low gain
709 name = "hPedlowgainDiff";
711 title = "Pedestals difference, low gain, module ";
713 fPedestalLowGainDiff.Add(new TProfile2D(name, title,
714 fColumns, 0.0, fColumns,
715 fRows, fRowMin, fRowMax,"s"));
717 //Pedestals, high gain
718 name = "hPedhighgainDiff";
720 title = "Pedestals difference, high gain, module ";
722 fPedestalHighGainDiff.Add(new TProfile2D(name, title,
723 fColumns, 0.0, fColumns,
724 fRows, fRowMin, fRowMax,"s"));
726 //LED Ref/Mon pedestals, low gain
727 name = "hPedestalLEDReflowgainDiff";
729 title = "Pedestal difference LEDRef, low gain, module ";
731 fPedestalLEDRefLowGainDiff.Add(new TProfile(name, title,
732 fLEDRefs, 0.0, fLEDRefs, "s"));
734 //LED Ref/Mon pedestals, high gain
735 name = "hPedestalLEDRefhighgainDiff";
737 title = "Pedestal difference LEDRef, high gain, module ";
739 fPedestalLEDRefHighGainDiff.Add(new TProfile(name, title,
740 fLEDRefs, 0.0, fLEDRefs, "s"));
742 //Peak-Pedestals, high gain
743 name = "hPeakMinusPedhighgainDiff";
745 title = "Peak-Pedestal difference, high gain, module ";
747 fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
748 fColumns, 0.0, fColumns,
749 fRows, fRowMin, fRowMax,"s"));
751 //Peak-Pedestals, low gain
752 name = "hPeakMinusPedlowgainDiff";
754 title = "Peak-Pedestal difference, low gain, module ";
756 fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
757 fColumns, 0.0, fColumns,
758 fRows, fRowMin, fRowMax,"s"));
760 //Pedestals, low gain
761 name = "hPedlowgainRatio";
763 title = "Pedestals ratio, low gain, module ";
765 fPedestalLowGainRatio.Add(new TProfile2D(name, title,
766 fColumns, 0.0, fColumns,
767 fRows, fRowMin, fRowMax,"s"));
769 //Pedestals, high gain
770 name = "hPedhighgainRatio";
772 title = "Pedestals ratio, high gain, module ";
774 fPedestalHighGainRatio.Add(new TProfile2D(name, title,
775 fColumns, 0.0, fColumns,
776 fRows, fRowMin, fRowMax,"s"));
778 //LED Ref/Mon pedestals, low gain
779 name = "hPedestalLEDReflowgainRatio";
781 title = "Pedestal ratio LEDRef, low gain, module ";
783 fPedestalLEDRefLowGainRatio.Add(new TProfile(name, title,
784 fLEDRefs, 0.0, fLEDRefs, "s"));
786 //LED Ref/Mon pedestals, high gain
787 name = "hPedestalLEDRefhighgainRatio";
789 title = "Pedestal ratio LEDRef, high gain, module ";
791 fPedestalLEDRefHighGainRatio.Add(new TProfile(name, title,
792 fLEDRefs, 0.0, fLEDRefs, "s"));
794 //Peak-Pedestals, low gain
795 name = "hPeakMinusPedlowgainRatio";
797 title = "Peak-Pedestal ratio, low gain, module ";
799 fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
800 fColumns, 0.0, fColumns,
801 fRows, fRowMin, fRowMax,"s"));
803 //Peak-Pedestals, high gain
804 name = "hPeakMinusPedhighgainRatio";
806 title = "Peak-Pedestal ratio, high gain, module ";
808 fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
809 fColumns, 0.0, fColumns,
810 fRows, fRowMin, fRowMax,"s"));
812 }//end for nModules create the histograms
815 //_____________________________________________________________________
816 void AliCaloCalibPedestal::ComputeDiffAndRatio()
818 // calculate differences and ratios relative to a reference
819 ValidateComparisonProfiles();//Make sure the comparison histos exist
822 return;//Return if the reference object isn't loaded
828 for (int i = 0; i < fModules; i++) {
829 //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
830 //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
831 for (int j = 0; j < fColumns; j++) {
832 for (int k = 0; k < fRows; k++) {
833 bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
835 if (fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) > 0) {
836 diff = GetPeakProfileHighGain(i)->GetBinContent(bin) - fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
837 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(bin, diff);
838 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
839 ratio = GetPeakProfileHighGain(i)->GetBinContent(bin) / fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
840 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinContent(bin, ratio);
841 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinEntries(bin, 1);
844 if (fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) > 0) {
845 diff = GetPeakProfileLowGain(i)->GetBinContent(bin) - fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
846 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(bin, diff);
847 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
848 ratio = GetPeakProfileLowGain(i)->GetBinContent(bin) / fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
849 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinContent(bin, ratio);
850 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinEntries(bin, 1);
853 if (fReference->GetPedProfileHighGain(i)->GetBinContent(bin) > 0) {
854 diff = GetPedProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
855 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(bin, diff);
856 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
857 ratio = GetPedProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
858 ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinContent(bin, ratio);
859 ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinEntries(bin, 1);
862 if (fReference->GetPedProfileLowGain(i)->GetBinContent(bin) > 0) {
863 diff = GetPedProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
864 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(bin, diff);
865 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
866 ratio = GetPedProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
867 ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinContent(bin, ratio);
868 ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinEntries(bin, 1);
874 // same for LED Ref/Mon channels
875 for (int j = 0; j <= fLEDRefs; j++) {
876 bin = j+1;//Note that we assume here that all histos have the same structure...
878 if (fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) > 0) {
879 diff = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
880 ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinContent(bin, diff);
881 ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinEntries(bin, 1);
882 ratio = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
883 ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinContent(bin, ratio);
884 ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinEntries(bin, 1);
887 if (fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) > 0) {
888 diff = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
889 ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinContent(bin, diff);
890 ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinEntries(bin, 1);
891 ratio = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
892 ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinContent(bin, ratio);
893 ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinEntries(bin, 1);
902 //_____________________________________________________________________
903 void AliCaloCalibPedestal::ComputeHotAndWarningTowers(const char * hotMapFile)
904 { // look for hot/noisy towers
906 char name[512];//Quite a long temp buffer, just in case the filename includes a path
909 snprintf(name, 512, "%s.txt", hotMapFile);
910 fout = new ofstream(name);
911 if (!fout->is_open()) {
913 fout = 0;//Set the pointer to empty if the file was not opened
917 for(int i = 0; i < fModules; i++){
919 //first we compute the peak-pedestal distribution for each supermodule...
920 if( GetPeakHighGainHisto(i)->GetEntries() > 0 ) {
921 double min = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMinimumBin());
922 double max = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMaximumBin());
923 TH1D *hPeakFit = new TH1D(Form("hFit_%d", i), Form("hFit_%d", i), (int)((max-min)*10), min-1, max+1);
925 for (int j = 1; j <= fColumns; j++) {
926 for (int k = 1; k <= fRows; k++) {
927 hPeakFit->Fill(GetPeakProfileHighGain(i)->GetBinContent(j, k));
931 //...and then we try to fit it
932 double mean = hPeakFit->GetMean();
933 double sigma = hPeakFit->GetRMS();
935 hPeakFit->Fit("gaus", "OQ", "", mean - 3*sigma, mean + 3*sigma);
936 mean = hPeakFit->GetFunction("gaus")->GetParameter(1);
937 sigma = hPeakFit->GetFunction("gaus")->GetParameter(2);
939 catch (const std::exception & e) {
940 printf("AliCaloCalibPedestal: TH1D PeakFit exception %s", e.what());
946 //Then we look for warm/hot towers
947 TH2F * hPeak2D = GetPeakHighGainHisto(i);
948 hPeak2D->GetYaxis()->SetRangeUser( fWarningThreshold, hPeak2D->GetYaxis()->GetBinUpEdge(hPeak2D->GetNbinsY()) );
952 for (int j = 1; j <= fColumns; j++) {
953 for (int k = 1; k <= fRows; k++) {
954 //we start looking for warm/warning towers...
955 // histogram x-axis index
956 idx = k-1 + fRows*(j-1); // this is what is used in the Fill call
957 hPeak2D->GetXaxis()->SetRangeUser(idx, idx);
958 warnCounter = (int) hPeak2D->Integral();
959 if(warnCounter > fNEvents * fWarningFraction) {
960 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kWarning);
961 /* printf("mod %d col %d row %d warnCounter %d - status %d\n",
962 i, j-1, k-1, warnCounter, (int) (kWarning)); */
964 //...then we look for hot ones (towers whose values are greater than mean + X*sigma)
965 if(GetPeakProfileHighGain(i)->GetBinContent(j, k) > mean + fHotSigma*sigma ) {
966 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kHot);
967 /* printf("mod %d col %d row %d binc %d - status %d\n",
968 i, j-1, k-1, (int)(GetPeakProfileHighGain(i)->GetBinContent(j, k)), (int) (kHot)); */
971 //Write the status to the hot/warm map file, if the file is open.
972 // module - column - row - status (1=dead, 2= warm/warning , 3 = hot, see .h file enum)
973 if (fout && ((TH2D*)fDeadMap[i])->GetBinContent(j, k) > 1) {
978 << ((TH2D*)fDeadMap[i])->GetBinContent(j, k) << endl; }
988 //_____________________________________________________________________
989 void AliCaloCalibPedestal::ComputeDeadTowers(const char * deadMapFile)
991 //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
997 char name[512];//Quite a long temp buffer, just in case the filename includes a path
1000 snprintf(name, 512, "%s.txt", deadMapFile);
1001 fout = new ofstream(name);
1002 snprintf(name, 512, "%sdiff.txt", deadMapFile);
1003 diff = new ofstream(name);
1004 if (!fout->is_open()) {
1006 fout = 0;//Set the pointer to empty if the file was not opened
1008 if (!diff->is_open()) {
1010 fout = 0;//Set the pointer to empty if the file was not opened
1014 for (int i = 0; i < fModules; i++) {
1015 if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
1016 for (int j = 1; j <= fColumns; j++) {
1017 for (int k = 1; k <= fRows; k++) {
1019 if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {//It's dead
1020 countTot++;//One more dead total
1026 << "0" << endl;//Write the status to the deadmap file, if the file is open.
1029 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= fDeadThreshold) {
1030 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
1031 countNew++;//This tower wasn't dead before!
1033 ( *diff) << i << " "
1037 << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
1041 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
1044 else { //It's ALIVE!!
1045 //Don't bother with writing the live ones.
1046 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {
1047 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
1048 countRes++; //This tower was dead before => it's a miracle! :P
1054 << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
1058 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
1063 }//end for j/columns
1064 }//end if GetEntries >= 0
1073 fDeadTowers = countTot;
1074 fNewDeadTowers = countNew;
1075 fResurrectedTowers = countRes;
1078 //_____________________________________________________________________
1079 Bool_t AliCaloCalibPedestal::IsBadChannel(int imod, int icol, int irow) const
1081 //Check if channel is dead or hot.
1082 Int_t status = (Int_t) ( ((TH2D*)fDeadMap[imod])->GetBinContent(icol,irow) );
1083 if(status == kAlive)
1090 //_____________________________________________________________________
1091 void AliCaloCalibPedestal::SetChannelStatus(int imod, int icol, int irow, int status)
1093 //Set status of channel dead, hot, alive ...
1094 ((TH2D*)fDeadMap[imod])->SetBinContent(icol, irow, status);