[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx

/**************************************************************************
 * 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: AliCaloCalibSignal.cxx $ */

//________________________________________________________________________
//
// 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
  fSignals = new AliCaloCalibSignal( fDetType );
  // AliCaloCalibSignal knows how many modules we have for PHOS or EMCAL
  fNumModules = fSignals->GetModules();
*/
// fed an event:
//  fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
// get some info:
//  fSignals->GetXXX..()
// etc.
//________________________________________________________________________

#include "TProfile.h"
#include "TFile.h"

#include "AliRawReader.h"
#include "AliRawEventHeaderBase.h"
#include "AliCaloRawStream.h"

//The include file
#include "AliCaloCalibSignal.h"

ClassImp(AliCaloCalibSignal)

using namespace std;

// variables for TTree filling; not sure if they should be static or not
static int fChannelNum;
static double fAmp;
static double fAvgAmp;
static double fRMS;

// ctor; initialize everything in order to avoid compiler warnings
// put some reasonable defaults
AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :  
  TObject(),
  fDetType(kNone),
  fColumns(0),
  fRows(0),
  fModules(0),
  fCaloString(),
  fMapping(NULL),
  fRunNumber(-1),
  fStartTime(0),
  fAmpCut(50),
  fReqFractionAboveAmpCutVal(0.8),
  fReqFractionAboveAmp(kTRUE),
  fHour(0),
  fLatestHour(0),
  fUseAverage(kTRUE),
  fSecInAverage(1800), 
  fNEvents(0),
  fNAcceptedEvents(0),
  fTreeAmpVsTime(NULL),
  fTreeAvgAmpVsTime(NULL)
{
  //Default constructor. First we set the detector-type related constants.
  if (detectorType == kPhos) {
    fColumns = fgkPhosCols;
    fRows = fgkPhosRows;
    fModules = fgkPhosModules;
    fCaloString = "PHOS";
  } 
  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;
    fCaloString = "EMCAL";
  }

  fDetType = detectorType;

  ResetInfo(); // trees and counters
} 

// dtor
//_____________________________________________________________________
AliCaloCalibSignal::~AliCaloCalibSignal()
{
  DeleteTrees();
}

//_____________________________________________________________________
void AliCaloCalibSignal::DeleteTrees()
{
  // delete what was created in the ctor (TTrees)
  if (fTreeAmpVsTime) delete fTreeAmpVsTime;
  if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
  // and reset pointers
  fTreeAmpVsTime = NULL;
  fTreeAvgAmpVsTime = NULL;

  return;
}

// copy ctor
//_____________________________________________________________________
AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
  TObject(sig),
  fDetType(sig.GetDetectorType()),
  fColumns(sig.GetColumns()),
  fRows(sig.GetRows()),
  fModules(sig.GetModules()),
  fCaloString(sig.GetCaloString()),
  fMapping(NULL), //! note that we are not copying the map info
  fRunNumber(sig.GetRunNumber()),
  fStartTime(sig.GetStartTime()),
  fAmpCut(sig.GetAmpCut()),
  fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
  fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
  fHour(sig.GetHour()),
  fLatestHour(sig.GetLatestHour()),
  fUseAverage(sig.GetUseAverage()),
  fSecInAverage(sig.GetSecInAverage()),
  fNEvents(sig.GetNEvents()),
  fNAcceptedEvents(sig.GetNAcceptedEvents()),
  fTreeAmpVsTime(NULL),
  fTreeAvgAmpVsTime(NULL)
{
  // also the TTree contents
  AddInfo(&sig);
}

// assignment operator; use copy ctor to make life easy..
//_____________________________________________________________________
AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
{
  // assignment operator; use copy ctor
  if (&source == this) return *this;

  new (this) AliCaloCalibSignal(source);
  return *this;
}

//_____________________________________________________________________
void AliCaloCalibSignal::CreateTrees()
{
  // initialize trees
  // first, regular version
  fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");

  fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
  fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
  fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");

  // then, average version
  fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");

  fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
  fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
  fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
  fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");

  return;
}

//_____________________________________________________________________
void AliCaloCalibSignal::ResetInfo()
{
  Zero(); // set all counters to 0
  DeleteTrees(); // delete previous stuff
  CreateTrees(); // and create some new ones
  return;
}

//_____________________________________________________________________
void AliCaloCalibSignal::Zero()
{
  // set all counters to 0; not cuts etc. though
  fHour = 0;
  fLatestHour = 0;
  fNEvents = 0;
  fNAcceptedEvents = 0;

  // Set the number of points for each Amp vs. Time graph to 0
  memset(fNHighGain, 0, sizeof(fNHighGain));
  memset(fNLowGain, 0, sizeof(fNLowGain));

  return;
}

//_____________________________________________________________________
Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
{
  int nAbove = 0;
    
  int TowerNum = 0;
  for (int i = 0; i<fModules; i++) {
    for (int j = 0; j<fColumns; j++) {
      for (int k = 0; k<fRows; k++) {
	TowerNum = GetTowerNum(i,j,k);
	if (AmpVal[TowerNum] > fAmpCut) { 
	  nAbove++;
	}
      }
    }
  }
  
  double fraction = (1.0*nAbove) / nTotChan;
  
  if (fraction > fReqFractionAboveAmpCutVal) {  
    return true;
  }
  else return false;
}

//_____________________________________________________________________
Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
{
  // add info from sig's TTrees to ours..
  TTree *sigAmp = sig->GetTreeAmpVsTime();
  TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();

  // we could try some merging via TList or what also as a more elegant approach
  // but I wanted with the stupid/simple and hopefully safe approach of looping
  // over what we want to add..

  // associate variables for sigAmp and sigAvgAmp:
  sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
  sigAmp->SetBranchAddress("fHour",&fHour);
  sigAmp->SetBranchAddress("fAmp",&fAmp);

  // loop over the trees.. note that since we use the same variables we should not need
  // to do any assignments between the getting and filling
  for (int i=0; i<sigAmp->GetEntries(); i++) {
    sigAmp->GetEntry(i);
    fTreeAmpVsTime->Fill();
  }

  sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
  sigAvgAmp->SetBranchAddress("fHour",&fHour);
  sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
  sigAvgAmp->SetBranchAddress("fRMS",&fRMS);

  for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
    sigAvgAmp->GetEntry(i);
    fTreeAvgAmpVsTime->Fill();
  }

  return kTRUE;//We hopefully succesfully added info from the supplied object
}

//_____________________________________________________________________
Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
{
  // if fMapping is NULL the rawstream will crate its own mapping
  AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);

  return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
}

//_____________________________________________________________________
Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream *in, AliRawEventHeaderBase *aliHeader)
{ 
  // 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

  // PHOS has more towers than EMCAL, so use PHOS numbers to set array sizes
  int AmpValHighGain[fgkMaxTowers];
  int AmpValLowGain[fgkMaxTowers];

  memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
  memset(AmpValLowGain, 0, sizeof(AmpValLowGain));

  int sample, isample = 0; //The sample temp, and the sample number in current event.
  int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal
  int gain = 0;
  
  // Number of Low and High gain channels for this event:
  int nLowChan = 0; 
  int nHighChan = 0; 

  int TowerNum = 0; // array index for TGraphs etc.

  // loop first to get the fraction of channels with amplitudes above cut
  while (in->Next()) {
    sample = in->GetSignal(); //Get the adc signal
    if (sample < min) min = sample;
    if (sample > max) max = sample;
    isample++;
    if ( isample >= in->GetTimeLength()) {
      //If we're here then we're done with this tower
      gain = 1 - in->IsLowGain();
      
      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("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos); 
      }

      // get tower number for AmpVal array
      TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 

      if (gain == 0) {
	// fill amplitude into the array	   
        AmpValLowGain[TowerNum] = max - min;
	nLowChan++;
      } 
      else if (gain==1) {//fill the high gain ones
	// fill amplitude into the array
	AmpValHighGain[TowerNum] = max - min;
	nHighChan++;
      }//end if gain

      
      max = fgkSampleMin; min = fgkSampleMax;
      isample = 0;
      
    }//End if end of tower
   
  }//end while, of stream
  
  // now check if it was a led event, only use high gain (that should be sufficient)
  if (fReqFractionAboveAmp) {
    bool ok = false;
    if (nHighChan > 0) { 
      ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan); 
    }
    if (!ok) return false; // skip event
  }

  fNAcceptedEvents++; // one more event accepted

  if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
    fStartTime = aliHeader->Get("Timestamp");
  }

  fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
  if (fLatestHour < fHour) {
    fLatestHour = fHour; 
  }
  
  // it is a led event, now fill TTree
  for(int i=0; i<fModules; i++){
    for(int j=0; j<fColumns; j++){
      for(int k=0; k<fRows; k++){
	
	TowerNum = GetTowerNum(i, j, k); 

	if(AmpValHighGain[TowerNum]) {
	  fAmp = AmpValHighGain[TowerNum];
	  fChannelNum = GetChannelNum(i,j,k,1);
	  fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
	  fNHighGain[TowerNum]++;
	}
	if(AmpValLowGain[TowerNum]) {
	  fAmp = AmpValLowGain[TowerNum];
	  fChannelNum = GetChannelNum(i,j,k,0);
	  fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
	  fNLowGain[TowerNum]++;
	}
      }
    }
  }
  
  return kTRUE;
}

//_____________________________________________________________________
Bool_t AliCaloCalibSignal::Save(TString fileName)
{
  //Saves all the TTrees to the designated file
  
  TFile destFile(fileName, "recreate");
  
  if (destFile.IsZombie()) {
    return kFALSE;
  }
  
  destFile.cd();

  // save the trees
  fTreeAmpVsTime->Write();
  if (fUseAverage) { 
    Analyze(); // get the latest and greatest averages
    fTreeAvgAmpVsTime->Write();
  }

  destFile.Close();
  
  return kTRUE;
}

//_____________________________________________________________________
Bool_t AliCaloCalibSignal::Analyze()
{
  // Fill the tree holding the average values
  if (!fUseAverage) { return kFALSE; }

  // Reset the average TTree if Analyze has already been called earlier,
  // meaning that the TTree could have been partially filled
  if (fTreeAvgAmpVsTime->GetEntries() > 0) {
    fTreeAvgAmpVsTime->Reset();
  }

  //0: setup variables for the TProfile plots that we'll use to do the averages
  int numProfBins = 0;
  double timeMin = 0;
  double timeMax = 0;
  if (fSecInAverage > 0) {
    numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
  }
  numProfBins += 2; // add extra buffer : first and last
  double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
  timeMin = - binSize;
  timeMax = timeMin + numProfBins*binSize;

  //1: set up TProfiles for the towers that had data
  TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
  memset(profile, 0, sizeof(profile));

  char name[200]; // for profile id and title
  int TowerNum = 0;

  for (int i = 0; i<fModules; i++) {
    for (int ic=0; ic<fColumns; ic++){
      for (int ir=0; ir<fRows; ir++) {

	TowerNum = GetTowerNum(i, ic, ir);
	// high gain
	if (fNHighGain[TowerNum] > 0) {
	  fChannelNum = GetChannelNum(i, ic, ir, 1); 
	  sprintf(name, "profileChan%d", fChannelNum);
	  profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
	}

	// same for low gain
	if (fNLowGain[TowerNum] > 0) {
	  fChannelNum = GetChannelNum(i, ic, ir, 0); 
	  sprintf(name, "profileChan%d", fChannelNum);
	  profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
	}

      } // rows
    } // columns
  } // modules

  //2: fill profiles by looping over tree
  // Set addresses for tree-readback also
  fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
  fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
  fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);

  for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
    fTreeAmpVsTime->GetEntry(ient);
    if (profile[fChannelNum]) { 
      // profile should always have been created above, for active channels
      profile[fChannelNum]->Fill(fHour, fAmp);
    }
  }

  // re-associating the branch addresses here seems to be needed for OK 'average' storage	    
  fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
  fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
  fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
  fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);

  //3: fill avg tree by looping over the profiles
  for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
    if (profile[fChannelNum]) { // profile was created
      if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
	for(int it=0; it<numProfBins; it++) {
	  if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
	    fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
	    fHour = profile[fChannelNum]->GetBinCenter(it+1);
	    fRMS = profile[fChannelNum]->GetBinError(it+1);
	    fTreeAvgAmpVsTime->Fill();
	  } // some entries for this bin
	} // loop over bins
      } // some entries for this profile
    } // profile exists  
  } // loop over all possible channels

  return kTRUE;
}
Commit	Line	Data
bd83f412	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	/* $Id: AliCaloCalibSignal.cxx $ */
	16
	17	//________________________________________________________________________
	18	//
	19	// A help class for monitoring and calibration tools: MOOD, AMORE etc.,
	20	// It can be created and used a la (ctor):
	21	/*
	22	//Create the object for making the histograms
	23	fSignals = new AliCaloCalibSignal( fDetType );
	24	// AliCaloCalibSignal knows how many modules we have for PHOS or EMCAL
	25	fNumModules = fSignals->GetModules();
	26	*/
	27	// fed an event:
	28	// fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
8bcca84c	29	// get some info:
8bcca84c	30	// fSignals->GetXXX..()
bd83f412	31	// etc.
	32	//________________________________________________________________________
	33
8bcca84c	34	#include "TProfile.h"
bd83f412	35	#include "TFile.h"
bd83f412	36
f4fc542c	37	#include "AliRawReader.h"
bd83f412	38	#include "AliRawEventHeaderBase.h"
	39	#include "AliCaloRawStream.h"
	40
	41	//The include file
	42	#include "AliCaloCalibSignal.h"
	43
	44	ClassImp(AliCaloCalibSignal)
	45
	46	using namespace std;
	47
8bcca84c	48	// variables for TTree filling; not sure if they should be static or not
	49	static int fChannelNum;
	50	static double fAmp;
	51	static double fAvgAmp;
	52	static double fRMS;
	53
bd83f412	54	// ctor; initialize everything in order to avoid compiler warnings
	55	// put some reasonable defaults
	56	AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
	57	TObject(),
	58	fDetType(kNone),
	59	fColumns(0),
	60	fRows(0),
	61	fModules(0),
f4fc542c	62	fCaloString(),
f4fc542c	63	fMapping(NULL),
bd83f412	64	fRunNumber(-1),
	65	fStartTime(0),
	66	fAmpCut(50),
	67	fReqFractionAboveAmpCutVal(0.8),
	68	fReqFractionAboveAmp(kTRUE),
	69	fHour(0),
	70	fLatestHour(0),
	71	fUseAverage(kTRUE),
	72	fSecInAverage(1800),
	73	fNEvents(0),
8bcca84c	74	fNAcceptedEvents(0),
	75	fTreeAmpVsTime(NULL),
	76	fTreeAvgAmpVsTime(NULL)
bd83f412	77	{
	78	//Default constructor. First we set the detector-type related constants.
	79	if (detectorType == kPhos) {
	80	fColumns = fgkPhosCols;
	81	fRows = fgkPhosRows;
	82	fModules = fgkPhosModules;
f4fc542c	83	fCaloString = "PHOS";
bd83f412	84	}
	85	else {
	86	//We'll just trust the enum to keep everything in line, so that if detectorType
	87	//isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
	88	//case, if someone intentionally gives another number
	89	fColumns = fgkEmCalCols;
	90	fRows = fgkEmCalRows;
	91	fModules = fgkEmCalModules;
f4fc542c	92	fCaloString = "EMCAL";
bd83f412	93	}
	94
	95	fDetType = detectorType;
	96
8bcca84c	97	ResetInfo(); // trees and counters
bd83f412	98	}
	99
	100	// dtor
	101	//_____________________________________________________________________
	102	AliCaloCalibSignal::~AliCaloCalibSignal()
	103	{
8bcca84c	104	DeleteTrees();
bd83f412	105	}
	106
	107	//_____________________________________________________________________
8bcca84c	108	void AliCaloCalibSignal::DeleteTrees()
bd83f412	109	{
8bcca84c	110	// delete what was created in the ctor (TTrees)
	111	if (fTreeAmpVsTime) delete fTreeAmpVsTime;
	112	if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
	113	// and reset pointers
	114	fTreeAmpVsTime = NULL;
	115	fTreeAvgAmpVsTime = NULL;
bd83f412	116
	117	return;
	118	}
	119
	120	// copy ctor
	121	//_____________________________________________________________________
	122	AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
	123	TObject(sig),
	124	fDetType(sig.GetDetectorType()),
	125	fColumns(sig.GetColumns()),
	126	fRows(sig.GetRows()),
	127	fModules(sig.GetModules()),
f4fc542c	128	fCaloString(sig.GetCaloString()),
f4fc542c	129	fMapping(NULL), //! note that we are not copying the map info
bd83f412	130	fRunNumber(sig.GetRunNumber()),
	131	fStartTime(sig.GetStartTime()),
	132	fAmpCut(sig.GetAmpCut()),
	133	fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
	134	fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
	135	fHour(sig.GetHour()),
	136	fLatestHour(sig.GetLatestHour()),
	137	fUseAverage(sig.GetUseAverage()),
	138	fSecInAverage(sig.GetSecInAverage()),
	139	fNEvents(sig.GetNEvents()),
8bcca84c	140	fNAcceptedEvents(sig.GetNAcceptedEvents()),
	141	fTreeAmpVsTime(NULL),
	142	fTreeAvgAmpVsTime(NULL)
bd83f412	143	{
8bcca84c	144	// also the TTree contents
bd83f412	145	AddInfo(&sig);
	146	}
	147
	148	// assignment operator; use copy ctor to make life easy..
	149	//_____________________________________________________________________
	150	AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
	151	{
	152	// assignment operator; use copy ctor
	153	if (&source == this) return *this;
	154
	155	new (this) AliCaloCalibSignal(source);
	156	return *this;
	157	}
	158
	159	//_____________________________________________________________________
8bcca84c	160	void AliCaloCalibSignal::CreateTrees()
bd83f412	161	{
8bcca84c	162	// initialize trees
	163	// first, regular version
	164	fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
	165
	166	fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
	167	fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
	168	fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
bd83f412	169
8bcca84c	170	// then, average version
	171	fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
	172
	173	fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
	174	fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
	175	fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
	176	fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
	177
	178	return;
bd83f412	179	}
	180
	181	//_____________________________________________________________________
8bcca84c	182	void AliCaloCalibSignal::ResetInfo()
bd83f412	183	{
bd83f412	184	Zero(); // set all counters to 0
8bcca84c	185	DeleteTrees(); // delete previous stuff
8bcca84c	186	CreateTrees(); // and create some new ones
bd83f412	187	return;
	188	}
	189
	190	//_____________________________________________________________________
	191	void AliCaloCalibSignal::Zero()
	192	{
8bcca84c	193	// set all counters to 0; not cuts etc. though
bd83f412	194	fHour = 0;
	195	fLatestHour = 0;
	196	fNEvents = 0;
	197	fNAcceptedEvents = 0;
8bcca84c	198
	199	// Set the number of points for each Amp vs. Time graph to 0
	200	memset(fNHighGain, 0, sizeof(fNHighGain));
	201	memset(fNLowGain, 0, sizeof(fNLowGain));
	202
bd83f412	203	return;
	204	}
	205
	206	//_____________________________________________________________________
	207	Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
	208	{
	209	int nAbove = 0;
	210
	211	int TowerNum = 0;
	212	for (int i = 0; i<fModules; i++) {
	213	for (int j = 0; j<fColumns; j++) {
	214	for (int k = 0; k<fRows; k++) {
	215	TowerNum = GetTowerNum(i,j,k);
	216	if (AmpVal[TowerNum] > fAmpCut) {
	217	nAbove++;
	218	}
	219	}
	220	}
	221	}
	222
	223	double fraction = (1.0*nAbove) / nTotChan;
	224
	225	if (fraction > fReqFractionAboveAmpCutVal) {
	226	return true;
	227	}
	228	else return false;
	229	}
	230
	231	//_____________________________________________________________________
	232	Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
	233	{
8bcca84c	234	// add info from sig's TTrees to ours..
	235	TTree *sigAmp = sig->GetTreeAmpVsTime();
	236	TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
	237
	238	// we could try some merging via TList or what also as a more elegant approach
	239	// but I wanted with the stupid/simple and hopefully safe approach of looping
	240	// over what we want to add..
	241
	242	// associate variables for sigAmp and sigAvgAmp:
	243	sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
	244	sigAmp->SetBranchAddress("fHour",&fHour);
	245	sigAmp->SetBranchAddress("fAmp",&fAmp);
	246
	247	// loop over the trees.. note that since we use the same variables we should not need
	248	// to do any assignments between the getting and filling
	249	for (int i=0; i<sigAmp->GetEntries(); i++) {
	250	sigAmp->GetEntry(i);
	251	fTreeAmpVsTime->Fill();
	252	}
4fb34e09	253
8bcca84c	254	sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
	255	sigAvgAmp->SetBranchAddress("fHour",&fHour);
	256	sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
	257	sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
bd83f412	258
8bcca84c	259	for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
	260	sigAvgAmp->GetEntry(i);
	261	fTreeAvgAmpVsTime->Fill();
	262	}
bd83f412	263
8bcca84c	264	return kTRUE;//We hopefully succesfully added info from the supplied object
bd83f412	265	}
bd83f412	266
f4fc542c	267	//_____________________________________________________________________
	268	Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
	269	{
	270	// if fMapping is NULL the rawstream will crate its own mapping
	271	AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
	272
	273	return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
	274	}
bd83f412	275
	276	//_____________________________________________________________________
	277	Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream in, AliRawEventHeaderBase aliHeader)
	278	{
	279	// Method to process=analyze one event in the data stream
	280	if (!in) return kFALSE; //Return right away if there's a null pointer
	281
	282	fNEvents++; // one more event
	283
	284	// PHOS has more towers than EMCAL, so use PHOS numbers to set array sizes
	285	int AmpValHighGain[fgkMaxTowers];
	286	int AmpValLowGain[fgkMaxTowers];
	287
	288	memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
	289	memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
	290
8cdb1ddb	291	int sample, isample = 0; //The sample temp, and the sample number in current event.
bd83f412	292	int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal
	293	int gain = 0;
	294
	295	// Number of Low and High gain channels for this event:
	296	int nLowChan = 0;
	297	int nHighChan = 0;
	298
	299	int TowerNum = 0; // array index for TGraphs etc.
	300
	301	// loop first to get the fraction of channels with amplitudes above cut
	302	while (in->Next()) {
	303	sample = in->GetSignal(); //Get the adc signal
	304	if (sample < min) min = sample;
	305	if (sample > max) max = sample;
8cdb1ddb	306	isample++;
8cdb1ddb	307	if ( isample >= in->GetTimeLength()) {
bd83f412	308	//If we're here then we're done with this tower
	309	gain = 1 - in->IsLowGain();
	310
	311	int arrayPos = in->GetModule(); //The modules are numbered starting from 0
	312	if (arrayPos >= fModules) {
	313	//TODO: return an error message, if appopriate (perhaps if debug>0?)
	314	return kFALSE;
	315	}
	316
	317	//Debug
	318	if (arrayPos < 0 \|\| arrayPos >= fModules) {
8bcca84c	319	printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
bd83f412	320	}
	321
	322	// get tower number for AmpVal array
	323	TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
	324
	325	if (gain == 0) {
	326	// fill amplitude into the array
	327	AmpValLowGain[TowerNum] = max - min;
	328	nLowChan++;
	329	}
	330	else if (gain==1) {//fill the high gain ones
	331	// fill amplitude into the array
	332	AmpValHighGain[TowerNum] = max - min;
	333	nHighChan++;
	334	}//end if gain
	335
	336
	337	max = fgkSampleMin; min = fgkSampleMax;
8cdb1ddb	338	isample = 0;
bd83f412	339
	340	}//End if end of tower
	341
	342	}//end while, of stream
	343
	344	// now check if it was a led event, only use high gain (that should be sufficient)
	345	if (fReqFractionAboveAmp) {
452729ca	346	bool ok = false;
	347	if (nHighChan > 0) {
	348	ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan);
	349	}
bd83f412	350	if (!ok) return false; // skip event
	351	}
	352
	353	fNAcceptedEvents++; // one more event accepted
	354
	355	if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
	356	fStartTime = aliHeader->Get("Timestamp");
	357	}
	358
	359	fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
	360	if (fLatestHour < fHour) {
	361	fLatestHour = fHour;
	362	}
	363
8bcca84c	364	// it is a led event, now fill TTree
	365	for(int i=0; i<fModules; i++){
	366	for(int j=0; j<fColumns; j++){
	367	for(int k=0; k<fRows; k++){
bd83f412	368
	369	TowerNum = GetTowerNum(i, j, k);
	370
	371	if(AmpValHighGain[TowerNum]) {
8bcca84c	372	fAmp = AmpValHighGain[TowerNum];
	373	fChannelNum = GetChannelNum(i,j,k,1);
	374	fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
	375	fNHighGain[TowerNum]++;
bd83f412	376	}
bd83f412	377	if(AmpValLowGain[TowerNum]) {
8bcca84c	378	fAmp = AmpValLowGain[TowerNum];
	379	fChannelNum = GetChannelNum(i,j,k,0);
	380	fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
	381	fNLowGain[TowerNum]++;
bd83f412	382	}
	383	}
	384	}
	385	}
	386
	387	return kTRUE;
	388	}
	389
	390	//_____________________________________________________________________
8bcca84c	391	Bool_t AliCaloCalibSignal::Save(TString fileName)
bd83f412	392	{
8bcca84c	393	//Saves all the TTrees to the designated file
bd83f412	394
	395	TFile destFile(fileName, "recreate");
	396
	397	if (destFile.IsZombie()) {
	398	return kFALSE;
	399	}
	400
	401	destFile.cd();
	402
8bcca84c	403	// save the trees
	404	fTreeAmpVsTime->Write();
	405	if (fUseAverage) {
	406	Analyze(); // get the latest and greatest averages
	407	fTreeAvgAmpVsTime->Write();
	408	}
	409
	410	destFile.Close();
	411
	412	return kTRUE;
	413	}
	414
	415	//_____________________________________________________________________
	416	Bool_t AliCaloCalibSignal::Analyze()
	417	{
	418	// Fill the tree holding the average values
	419	if (!fUseAverage) { return kFALSE; }
	420
	421	// Reset the average TTree if Analyze has already been called earlier,
	422	// meaning that the TTree could have been partially filled
	423	if (fTreeAvgAmpVsTime->GetEntries() > 0) {
	424	fTreeAvgAmpVsTime->Reset();
	425	}
	426
	427	//0: setup variables for the TProfile plots that we'll use to do the averages
bd83f412	428	int numProfBins = 0;
	429	double timeMin = 0;
	430	double timeMax = 0;
8bcca84c	431	if (fSecInAverage > 0) {
	432	numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
	433	}
	434	numProfBins += 2; // add extra buffer : first and last
	435	double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
	436	timeMin = - binSize;
	437	timeMax = timeMin + numProfBins*binSize;
	438
	439	//1: set up TProfiles for the towers that had data
	440	TProfile * profile[fgkMaxTowers2]; // 2 is since we include both high and low gains
	441	memset(profile, 0, sizeof(profile));
	442
	443	char name[200]; // for profile id and title
	444	int TowerNum = 0;
	445
	446	for (int i = 0; i<fModules; i++) {
	447	for (int ic=0; ic<fColumns; ic++){
	448	for (int ir=0; ir<fRows; ir++) {
	449
	450	TowerNum = GetTowerNum(i, ic, ir);
	451	// high gain
	452	if (fNHighGain[TowerNum] > 0) {
	453	fChannelNum = GetChannelNum(i, ic, ir, 1);
	454	sprintf(name, "profileChan%d", fChannelNum);
	455	profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
	456	}
	457
	458	// same for low gain
	459	if (fNLowGain[TowerNum] > 0) {
	460	fChannelNum = GetChannelNum(i, ic, ir, 0);
	461	sprintf(name, "profileChan%d", fChannelNum);
	462	profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
	463	}
	464
	465	} // rows
	466	} // columns
	467	} // modules
	468
	469	//2: fill profiles by looping over tree
	470	// Set addresses for tree-readback also
	471	fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
	472	fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
	473	fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
	474
	475	for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
	476	fTreeAmpVsTime->GetEntry(ient);
	477	if (profile[fChannelNum]) {
	478	// profile should always have been created above, for active channels
	479	profile[fChannelNum]->Fill(fHour, fAmp);
bd83f412	480	}
bd83f412	481	}
bd83f412	482
8bcca84c	483	// re-associating the branch addresses here seems to be needed for OK 'average' storage
	484	fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
	485	fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
	486	fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
	487	fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
	488
	489	//3: fill avg tree by looping over the profiles
	490	for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
	491	if (profile[fChannelNum]) { // profile was created
	492	if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
	493	for(int it=0; it<numProfBins; it++) {
	494	if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
	495	fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
	496	fHour = profile[fChannelNum]->GetBinCenter(it+1);
	497	fRMS = profile[fChannelNum]->GetBinError(it+1);
	498	fTreeAvgAmpVsTime->Fill();
	499	} // some entries for this bin
	500	} // loop over bins
	501	} // some entries for this profile
	502	} // profile exists
	503	} // loop over all possible channels
	504
bd83f412	505	return kTRUE;
bd83f412	506	}