+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+/* $Id$ */
+
+//________________________________________________________________________
+//
+// A help class for monitoring and calibration tools: MOOD, AMORE etc.,
+// It can be created and used a la (ctor):
+/*
+ //Create the object for making the histograms
+ fPedestals = new AliCaloCalibPedestal( fDetType );
+ // AliCaloCalibPedestal knows how many modules we have for PHOS or EMCAL
+ fNumModules = fPedestals->GetModules();
+*/
+// fed an event:
+// fPedestals->ProcessEvent(fCaloRawStream);
+// asked to draw histograms:
+// fPedestals->GetDeadMap(i)->Draw("col");
+// or
+// fPedestals->GetPeakProfileHighGainRatio((i < fNumModules) ? i : fVisibleModule)->Draw("colz");
+// etc.
+// The pseudo-code examples above were from the first implementation in MOOD (summer 2007).
+//________________________________________________________________________
+
#include "TH1.h"
#include "TFile.h"
-
#include <fstream>
+#include <sstream>
#include <iostream>
+#include <cmath>
+
+#include "AliCaloRawStreamV3.h"
//The include file
#include "AliCaloCalibPedestal.h"
TObject(),
fPedestalLowGain(),
fPedestalHighGain(),
+ fPedestalRMSLowGain(),
+ fPedestalRMSHighGain(),
fPeakMinusPedLowGain(),
fPeakMinusPedHighGain(),
fPedestalLowGainDiff(),
fPeakMinusPedLowGainRatio(),
fPeakMinusPedHighGainRatio(),
fDeadMap(),
+ fNEvents(0),
+ fNChanFills(0),
fDeadTowers(0),
fNewDeadTowers(0),
fResurrectedTowers(0),
fColumns(0),
fRows(0),
fModules(0),
- fRunNumber(-1)
+ fRowMin(0),
+ fRowMax(0),
+ fRowMultiplier(0),
+ fCaloString(),
+ fMapping(NULL),
+ fRunNumber(-1),
+ fSelectPedestalSamples(kTRUE),
+ fFirstPedestalSample(0),
+ fLastPedestalSample(15)
{
//Default constructor. First we set the detector-type related constants.
if (detectorType == kPhos) {
fColumns = fgkPhosCols;
fRows = fgkPhosRows;
fModules = fgkPhosModules;
+ fCaloString = "PHOS";
+ fRowMin = -1*fRows;
+ fRowMax = 0;
+ fRowMultiplier = -1;
}
else {
//We'll just trust the enum to keep everything in line, so that if detectorType
//isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
//case, if someone intentionally gives another number
- fColumns = fgkEmCalCols;
- fRows = fgkEmCalRows;
- fModules = fgkEmCalModules;
+ fColumns = AliEMCALGeoParams::fgkEMCALCols;
+ fRows = AliEMCALGeoParams::fgkEMCALRows;
+ fModules = AliEMCALGeoParams::fgkEMCALModules;
+ fCaloString = "EMCAL";
+ fRowMin = 0;
+ fRowMax = fRows;
+ fRowMultiplier = 1;
}
fDetType = detectorType;
title += i;
fPedestalLowGain.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Pedestals, high gain
name = "hPedhighgain";
title += i;
fPedestalHighGain.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
+ //All Samples, low gain
+ name = "hPedestalRMSlowgain";
+ name += i;
+ title = "Pedestal RMS, low gain, module ";
+ title += i;
+ fPedestalRMSLowGain.Add(new TProfile2D(name, title,
+ fColumns, 0.0, fColumns,
+ fRows, fRowMin, fRowMax,"s"));
+
+ //All Samples, high gain
+ name = "hPedestalRMShighgain";
+ name += i;
+ title = "Pedestal RMS, high gain, module ";
+ title += i;
+ fPedestalRMSHighGain.Add(new TProfile2D(name, title,
+ fColumns, 0.0, fColumns,
+ fRows, fRowMin, fRowMax,"s"));
//Peak-Pedestals, low gain
name = "hPeakMinusPedlowgain";
title += i;
fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Peak-Pedestals, high gain
name = "hPeakMinusPedhighgain";
title += i;
fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
name = "hDeadMap";
name += i;
title = "Dead map, module ";
title += i;
fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax));
}//end for nModules create the histograms
//Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
fPedestalLowGain.Compress();
fPedestalHighGain.Compress();
+ fPedestalRMSLowGain.Compress();
+ fPedestalRMSHighGain.Compress();
fPeakMinusPedLowGain.Compress();
fPeakMinusPedHighGain.Compress();
fDeadMap.Compress();
TObject(ped),
fPedestalLowGain(),
fPedestalHighGain(),
+ fPedestalRMSLowGain(),
+ fPedestalRMSHighGain(),
fPeakMinusPedLowGain(),
fPeakMinusPedHighGain(),
fPedestalLowGainDiff(),
fPeakMinusPedLowGainRatio(),
fPeakMinusPedHighGainRatio(),
fDeadMap(),
+ fNEvents(ped.GetNEvents()),
+ fNChanFills(ped.GetNChanFills()),
fDeadTowers(ped.GetDeadTowerCount()),
fNewDeadTowers(ped.GetDeadTowerNew()),
fResurrectedTowers(ped.GetDeadTowerResurrected()),
fColumns(ped.GetColumns()),
fRows(ped.GetRows()),
fModules(ped.GetModules()),
- fRunNumber(ped.GetRunNumber())
+ fRowMin(ped.GetRowMin()),
+ fRowMax(ped.GetRowMax()),
+ fRowMultiplier(ped.GetRowMultiplier()),
+ fCaloString(ped.GetCaloString()),
+ fMapping(NULL), //! note that we are not copying the map info
+ fRunNumber(ped.GetRunNumber()),
+ fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
+ fFirstPedestalSample(ped.GetFirstPedestalSample()),
+ fLastPedestalSample(ped.GetLastPedestalSample())
{
// Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
//DS: this has not really been tested yet..
for (int i = 0; i < fModules; i++) {
fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
+ fPedestalRMSLowGain.Add( ped.GetPedRMSProfileLowGain(i) );
+ fPedestalRMSHighGain.Add( ped.GetPedRMSProfileHighGain(i) );
fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
//Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
fPedestalLowGain.Compress();
fPedestalHighGain.Compress();
+ fPedestalRMSLowGain.Compress();
+ fPedestalRMSHighGain.Compress();
fPeakMinusPedLowGain.Compress();
fPeakMinusPedHighGain.Compress();
fDeadMap.Compress();
//_____________________________________________________________________
AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
{
+ // assignment operator; use copy ctor
if (&source == this) return *this;
new (this) AliCaloCalibPedestal(source);
//_____________________________________________________________________
void AliCaloCalibPedestal::Reset()
{
+ // Reset all arrays/histograms
for (int i = 0; i < fModules; i++) {
GetPedProfileLowGain(i)->Reset();
GetPedProfileHighGain(i)->Reset();
GetPeakProfileHighGainRatio(i)->Reset();
}
}
+ fNEvents = 0;
+ fNChanFills = 0;
fDeadTowers = 0;
fNewDeadTowers = 0;
fResurrectedTowers = 0;
//To think about: should fReference be deleted too?... let's not do it this time, at least...
}
+// Parameter/cut handling
+//_____________________________________________________________________
+void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile)
+{
+ static const string delimitor("::");
+
+ // open, check input file
+ ifstream in( parameterFile );
+ if( !in ) {
+ printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
+ return;
+ }
+
+ // Note: this method is a bit more complicated than it really has to be
+ // - allowing for multiple entries per line, arbitrary order of the
+ // different variables etc. But I wanted to try and do this in as
+ // correct a C++ way as I could (as an exercise).
+
+ // read in
+ char readline[1024];
+ while ((in.rdstate() & ios::failbit) == 0 ) {
+
+ // Read into the raw char array and then construct a string
+ // to do the searching
+ in.getline(readline, 1024);
+ istringstream s(readline);
+
+ while ( ( s.rdstate() & ios::failbit ) == 0 ) {
+
+ string keyValue;
+ s >> keyValue;
+
+ // check stream status
+ if( s.rdstate() & ios::failbit ) break;
+
+ // skip rest of line if comments found
+ if( keyValue.substr( 0, 2 ) == "//" ) break;
+
+ // look for "::" in keyValue pair
+ size_t position = keyValue.find( delimitor );
+ if( position == string::npos ) {
+ printf("wrong format for key::value pair: %s\n", keyValue.c_str());
+ }
+
+ // split keyValue pair
+ string key( keyValue.substr( 0, position ) );
+ string value( keyValue.substr( position+delimitor.size(),
+ keyValue.size()-delimitor.size() ) );
+
+ // check value does not contain a new delimitor
+ if( value.find( delimitor ) != string::npos ) {
+ printf("wrong format for key::value pair: %s\n", keyValue.c_str());
+ }
+
+ // debug: check key value pair
+ // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
+
+ // if the key matches with something we expect, we assign the new value
+ istringstream iss(value);
+ // the comparison strings defined at the beginning of this method
+ if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") ) {
+ printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
+
+ if (key == "fFirstPedestalSample") {
+ iss >> fFirstPedestalSample;
+ }
+ else if (key == "fLastPedestalSample") {
+ iss >> fLastPedestalSample;
+ }
+ } // some match
+
+ }
+ }
+
+ in.close();
+ return;
+
+}
+
+//_____________________________________________________________________
+void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
+{
+ static const string delimitor("::");
+ ofstream out( parameterFile );
+ out << "// " << parameterFile << endl;
+ out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
+ out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
+
+ out.close();
+ return;
+}
+
+//_____________________________________________________________________
+Bool_t AliCaloCalibPedestal::AddInfo(const AliCaloCalibPedestal *ped)
+{
+ // just do this for the basic histograms/profiles that get filled in ProcessEvent
+ // may not have data for all modules, but let's just Add everything..
+ for (int i = 0; i < fModules; i++) {
+ GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
+ GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
+ GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
+ GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
+ }//end for nModules
+
+ // DeadMap; Diff profiles etc would need to be redone after this operation
+
+ return kTRUE;//We succesfully added info from the supplied object
+}
+
+//_____________________________________________________________________
+Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
+{
+ // if fMapping is NULL the rawstream will crate its own mapping
+ AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
+ return ProcessEvent(&rawStream);
+}
+
//_____________________________________________________________________
-Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
+Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
{
+ // Method to process=analyze one event in the data stream
if (!in) return kFALSE; //Return right away if there's a null pointer
+ fNEvents++; // one more event
- in->SetOldRCUFormat(kTRUE);
-
- int sample, i = 0; //The sample temp, and the sample number in current event.
- int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
+ // indices for the reading
+ int sample = 0;
int gain = 0;
-
- while (in->Next()) {
- sample = in->GetSignal(); //Get the adc signal
- if (sample < min) min = sample;
- if (sample > max) max = sample;
- i++;
- if ( i >= in->GetTimeLength()) {
- //If we're here then we're done with this tower
- gain = 1 - in->IsLowGain();
+ int time = 0;
+ int i = 0; // sample counter
+ int startBin = 0;
+
+ // start loop over input stream
+ while (in->NextDDL()) {
+ while (in->NextChannel()) {
+
+ // counters
+ int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
+ int nsamples = 0;
+
+ // for the pedestal calculation
+ int sampleSum = 0; // sum of samples
+ int squaredSampleSum = 0; // sum of samples squared
+ int nSum = 0; // number of samples in sum
+ // calc. quantities
+ double mean = 0, squaredMean = 0, rms = 0;
+
+ while (in->NextBunch()) {
+ const UShort_t *sig = in->GetSignals();
+ startBin = in->GetStartTimeBin();
+ nsamples += in->GetBunchLength();
+ for (i = 0; i < in->GetBunchLength(); i++) {
+ sample = sig[i];
+ time = startBin--;
+
+ // check if it's a min or max value
+ if (sample < min) min = sample;
+ if (sample > max) max = sample;
+
+ // should we add it for the pedestal calculation?
+ if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
+ !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
+ sampleSum += sample;
+ squaredSampleSum += sample*sample;
+ nSum++;
+ }
+
+ } // loop over samples in bunch
+ } // loop over bunches
+
+ if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
+
+ // calculate pedesstal estimate: mean of possibly selected samples
+ if (nSum > 0) {
+ mean = sampleSum / (1.0 * nSum);
+ squaredMean = squaredSampleSum / (1.0 * nSum);
+ // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
+ rms = sqrt(squaredMean - mean*mean);
+ }
+ else {
+ mean = 0;
+ squaredMean = 0;
+ rms = 0;
+ }
+
+ // we're done with the calc. for this channel; let's prepare to fill histo
+ gain = -1; // init to not valid value
+ if ( in->IsLowGain() ) {
+ gain = 0;
+ }
+ else if ( in->IsHighGain() ) {
+ gain = 1;
+ }
+ // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
int arrayPos = in->GetModule(); //The modules are numbered starting from 0
if (arrayPos >= fModules) {
//TODO: return an error message, if appopriate (perhaps if debug>0?)
return kFALSE;
- }
-
+ }
//Debug
if (arrayPos < 0 || arrayPos >= fModules) {
printf("Oh no: arrayPos = %i.\n", arrayPos);
}
-
+
+ fNChanFills++; // one more channel found, and profile to be filled
//NOTE: coordinates are (column, row) for the profiles
if (gain == 0) {
//fill the low gain histograms
- ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
- ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
+ ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
+ if (nSum>0) { // only fill pedestal info in case it could be calculated
+ ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
+ ((TProfile2D*)fPedestalRMSLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
+ }
}
- else {//fill the high gain ones
- ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
- ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
- }//end if gain
-
- max = fgkSampleMin; min = fgkSampleMax;
- i = 0;
-
- }//End if end of tower
-
- }//end while, of stream
-
+ else if (gain == 1) {
+ //fill the high gain ones
+ ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
+ if (nSum>0) { // only fill pedestal info in case it could be calculated
+ ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
+ ((TProfile2D*)fPedestalRMSHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
+ }
+ }//end if valid gain
+
+ } // nsamples>0 check, some data found for this channel; not only trailer/header
+ }// end while over channel
+ }//end while over DDL's, of input stream
+
+ in->Reset(); // just in case the next customer forgets to check if the stream was reset..
+
return kTRUE;
}
}
if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
fPedestalHighGain[i]->Write();
+ Printf("save %d", i);
+ }
+ if( ((TProfile2D *)fPedestalRMSLowGain[i])->GetEntries() || saveEmptyHistos) {
+ fPedestalRMSLowGain[i]->Write();
+ }
+ if( ((TProfile2D *)fPedestalRMSHighGain[i])->GetEntries() || saveEmptyHistos) {
+ fPedestalRMSHighGain[i]->Write();
}
}
//_____________________________________________________________________
void AliCaloCalibPedestal::ValidateComparisonProfiles()
{
+ //Make sure the comparison histos exist
if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
//the same time
title += i;
fPedestalLowGainDiff.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Pedestals, high gain
name = "hPedhighgainDiff";
title += i;
fPedestalHighGainDiff.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
-
- //Peak-Pedestals, low gain
- name = "hPeakMinusPedlowgainDiff";
- name += i;
- title = "Peak-Pedestal difference, low gain, module ";
- title += i;
- fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
- fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
-
+ fRows, fRowMin, fRowMax,"s"));
+
//Peak-Pedestals, high gain
name = "hPeakMinusPedhighgainDiff";
name += i;
title += i;
fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Pedestals, low gain
name = "hPedlowgainRatio";
title += i;
fPedestalLowGainRatio.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Pedestals, high gain
name = "hPedhighgainRatio";
title += i;
fPedestalHighGainRatio.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Peak-Pedestals, low gain
name = "hPeakMinusPedlowgainRatio";
title += i;
fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
//Peak-Pedestals, high gain
name = "hPeakMinusPedhighgainRatio";
title += i;
fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
fColumns, 0.0, fColumns,
- fRows, -fRows, 0.0));
+ fRows, fRowMin, fRowMax,"s"));
}//end for nModules create the histograms
}
//_____________________________________________________________________
void AliCaloCalibPedestal::ComputeDiffAndRatio()
{
-
+ // calculate differences and ratios relative to a reference
ValidateComparisonProfiles();//Make sure the comparison histos exist
if (!fReference) {
//_____________________________________________________________________
void AliCaloCalibPedestal::ComputeDeadTowers(int threshold, const char * deadMapFile)
-{//Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
+{
+ //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
int countTot = 0;
int countNew = 0;
int countRes = 0;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff