[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);
// asked to draw graphs or profiles:
//  fSignals->GetGraphAmpVsTimeHighGain(module,column,row)->Draw("ap");
// or
//  fSignals->GetProfAmpVsTimeHighGain(module,column,row)->Draw();
// etc.
//________________________________________________________________________

#include "TFile.h"

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

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

ClassImp(AliCaloCalibSignal)

using namespace std;

// 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)
{
  //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;

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

  CreateGraphs(); // set up the TGraphs

  // init TProfiles to NULL=0 also
  memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
  memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
} 

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

//_____________________________________________________________________
void AliCaloCalibSignal::ClearObjects()
{
  // delete what was created in the ctor (TGraphs), and possible later (TProfiles)
  for (int i=0; i<fgkMaxTowers; i++) {
    if ( fGraphAmpVsTimeHighGain[i] ) { delete fGraphAmpVsTimeHighGain[i]; }
    if ( fGraphAmpVsTimeLowGain[i] ) { delete fGraphAmpVsTimeLowGain[i]; }
    if ( fProfAmpVsTimeHighGain[i] ) { delete fProfAmpVsTimeHighGain[i]; }
    if ( fProfAmpVsTimeLowGain[i] ) { delete fProfAmpVsTimeLowGain[i]; }
  }
  // set pointers
  memset(fGraphAmpVsTimeHighGain, 0, sizeof(fGraphAmpVsTimeHighGain));
  memset(fGraphAmpVsTimeLowGain, 0, sizeof(fGraphAmpVsTimeLowGain));
  memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
  memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));

  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())
{
  // also the TGraph 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::CreateGraphs()
{
  //Then, loop for the requested number of modules
  TString title, name;
  for (int i = 0; i < fModules; i++) {
    
    // Amplitude vs. Time graph for each channel
    for(int ic=0;ic < fColumns;ic++){
      for(int ir=0;ir < fRows;ir++){
	
	int id = GetTowerNum(i, ic, ir);
	
	// high gain graph
	name = "fGraphAmpVsTimeHighGain_"; name += i;
	name += "_"; name += ic;
	name += "_"; name += ir;
	title = "Amp vs. Time High Gain Mod "; title += i;
	title += " Col "; title += ic;
	title += " Row "; title += ir;
	
	fGraphAmpVsTimeHighGain[id] = new TGraph();
	fGraphAmpVsTimeHighGain[id]->SetName(name);
	fGraphAmpVsTimeHighGain[id]->SetTitle(title);
	fGraphAmpVsTimeHighGain[id]->SetMarkerStyle(20);
	
	// Low Gain
	name = "fGraphAmpVsTimeLowGain_"; name += i;
	name += "_"; name += ic;
	name += "_"; name += ir;
	title = "Amp vs. Time Low Gain Mod "; title += i;
	title += " Col "; title += ic;
	title += " Row "; title += ir;
	
	fGraphAmpVsTimeLowGain[id] = new TGraph();
	fGraphAmpVsTimeLowGain[id]->SetName(name);
	fGraphAmpVsTimeLowGain[id]->SetTitle(title);
	fGraphAmpVsTimeLowGain[id]->SetMarkerStyle(20);
	
      }
    }

  }//end for nModules 
}

//_____________________________________________________________________
void AliCaloCalibSignal::Reset()
{
  Zero(); // set all counters to 0
  ClearObjects(); // delete previous TGraphs and TProfiles
  CreateGraphs(); // 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;
  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)
{
  // just do this for the basic graphs/profiles that get filled in ProcessEvent
  // may not have data for all channels, but let's just Add everything..
  // Note: this method will run into problems with TProfile adding if the binning of
  // the local profiles is not the same as those provided by the argument *sig..
  int numGraphPoints = 0;
  int id = 0;
  int ip = 0;
  for (int i = 0; i < fModules; i++) {
    for (int j = 0; j < fColumns; j++) {
      for (int k = 0; k < fRows; k++) {
	
	id = GetTowerNum(i,j,k);

	if(fUseAverage){ // add to Profiles
 	  if (sig->GetProfAmpVsTimeHighGain(id)) {
	    GetProfAmpVsTimeHighGain(id)->Add(sig->GetProfAmpVsTimeHighGain(id));
	  }
 	  if (sig->GetProfAmpVsTimeLowGain(id)) {
	    GetProfAmpVsTimeLowGain(id)->Add(sig->GetProfAmpVsTimeLowGain(id));
	  }
	}
	else{ // add to Graphs	  
	  // high gain
	  numGraphPoints= sig->GetGraphAmpVsTimeHighGain(id)->GetN();
	  if (numGraphPoints > 0) {
	    // get the values
	    double *graphX = sig->GetGraphAmpVsTimeHighGain(id)->GetX();
	    double *graphY = sig->GetGraphAmpVsTimeHighGain(id)->GetY();
	    for(ip=0; ip < numGraphPoints; ip++){
	      fGraphAmpVsTimeHighGain[id]->SetPoint(fNHighGain[id]++,graphX[ip],graphY[ip]);
	    }
	  }
	  // low gain
	  numGraphPoints= sig->GetGraphAmpVsTimeLowGain(id)->GetN();
	  if (numGraphPoints > 0) {
	    // get the values
	    double *graphX = sig->GetGraphAmpVsTimeLowGain(id)->GetX();
	    double *graphY = sig->GetGraphAmpVsTimeLowGain(id)->GetY();
	    for(ip=0; ip < numGraphPoints; ip++){
	      fGraphAmpVsTimeLowGain[id]->SetPoint(fNLowGain[id]++,graphX[ip],graphY[ip]);
	    }
	  }

	}

      }//end for nModules 
    }//end for nColumns
  }//end for nRows

  return kTRUE;//We 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("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 graphs (maybe profiles later)
  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]) {
	  fGraphAmpVsTimeHighGain[TowerNum]->SetPoint(fNHighGain[TowerNum]++,fHour,AmpValHighGain[TowerNum]);
	}
	if(AmpValLowGain[TowerNum]) {
 	  fGraphAmpVsTimeLowGain[TowerNum]->SetPoint(fNLowGain[TowerNum]++,fHour,AmpValLowGain[TowerNum]);
	}
      }
    }
  }
  
  return kTRUE;
}

//_____________________________________________________________________
void AliCaloCalibSignal::CreateProfile(int imod, int ic, int ir, int towerId, int gain,
				       int nbins, double min, double max)
{ //! create/setup a TProfile
  TString title, name;   
  if (gain == 0) { 
    name = "fProfAmpVsTimeLowGain_";   
    title = "Amp vs. Time Low Gain Mod "; 
  } 
  else if (gain == 1) { 
    name = "fProfAmpVsTimeHighGain_"; 
    title = "Amp vs. Time High Gain Mod "; 
  } 
  name += imod;
  name += "_"; name += ic;
  name += "_"; name += ir;
  title += imod;
  title += " Col "; title += ic;
  title += " Row "; title += ir;
	    
  // use "s" option for RMS
  if (gain==0) { 
    fProfAmpVsTimeLowGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
  }
  else if (gain==1) {
    fProfAmpVsTimeHighGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
  }

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

  // setup variables for the TProfile plot
  int numProfBins = 0;
  double timeMin = 0;
  double timeMax = 0;
  if (fUseAverage) {
    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;
  }

  int numGraphPoints= 0;
  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); 

	// 1st: high gain
	numGraphPoints= fGraphAmpVsTimeHighGain[TowerNum]->GetN();
	if( numGraphPoints>0 || saveEmptyGraphs) {
	  
	  // average the graphs points over time if requested and put them in a profile plot
	  if(fUseAverage && numGraphPoints>0) {
	    
	    // get the values
	    double *graphX = fGraphAmpVsTimeHighGain[TowerNum]->GetX();
	    double *graphY = fGraphAmpVsTimeHighGain[TowerNum]->GetY();

	    // create the TProfile: 1 is for High gain	    	    
	    CreateProfile(i, ic, ir, TowerNum, 1,
			  numProfBins, timeMin, timeMax);

	    // loop over graph points and fill profile
	    for(int ip=0; ip < numGraphPoints; ip++){
	      fProfAmpVsTimeHighGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
	    }
	    
	    fProfAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
	    fProfAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
	    fProfAmpVsTimeHighGain[TowerNum]->Write();

	  }
	   else{
	     //otherwise, just save the graphs and forget the profiling
	     fGraphAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
	     fGraphAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
	     fGraphAmpVsTimeHighGain[TowerNum]->Write();
	   }
	  
	} // low gain graph info should be saved in one form or another
	
	// 2nd: now go to the low gain case
	numGraphPoints= fGraphAmpVsTimeLowGain[TowerNum]->GetN();
	if( numGraphPoints>0 || saveEmptyGraphs) {
	  
	  // average the graphs points over time if requested and put them in a profile plot
	  if(fUseAverage && numGraphPoints>0) {
	    
	    double *graphX = fGraphAmpVsTimeLowGain[TowerNum]->GetX();
	    double *graphY = fGraphAmpVsTimeLowGain[TowerNum]->GetY();
	    
	    // create the TProfile: 0 is for Low gain	    
	    CreateProfile(i, ic, ir, TowerNum, 0,
			  numProfBins, timeMin, timeMax);

	    // loop over graph points and fill profile
	    for(int ip=0; ip < numGraphPoints; ip++){
	      fProfAmpVsTimeLowGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
	    }
	    
	    fProfAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
	    fProfAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
	    fProfAmpVsTimeLowGain[TowerNum]->Write();

	  }
	  else{
	     //otherwise, just save the graphs and forget the profiling
	    fGraphAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
	    fGraphAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
	    fGraphAmpVsTimeLowGain[TowerNum]->Write();
	  }
	  
	} // low gain graph info should be saved in one form or another

      } // fRows
    } // fColumns

  } // fModules
  destFile.Close();
  
  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);
	29	// asked to draw graphs or profiles:
	30	// fSignals->GetGraphAmpVsTimeHighGain(module,column,row)->Draw("ap");
	31	// or
	32	// fSignals->GetProfAmpVsTimeHighGain(module,column,row)->Draw();
	33	// etc.
	34	//________________________________________________________________________
	35
	36	#include "TFile.h"
	37
f4fc542c	38	#include "AliRawReader.h"
bd83f412	39	#include "AliRawEventHeaderBase.h"
	40	#include "AliCaloRawStream.h"
	41
	42	//The include file
	43	#include "AliCaloCalibSignal.h"
	44
	45	ClassImp(AliCaloCalibSignal)
	46
	47	using namespace std;
	48
	49	// ctor; initialize everything in order to avoid compiler warnings
	50	// put some reasonable defaults
	51	AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
	52	TObject(),
	53	fDetType(kNone),
	54	fColumns(0),
	55	fRows(0),
	56	fModules(0),
f4fc542c	57	fCaloString(),
f4fc542c	58	fMapping(NULL),
bd83f412	59	fRunNumber(-1),
	60	fStartTime(0),
	61	fAmpCut(50),
	62	fReqFractionAboveAmpCutVal(0.8),
	63	fReqFractionAboveAmp(kTRUE),
	64	fHour(0),
	65	fLatestHour(0),
	66	fUseAverage(kTRUE),
	67	fSecInAverage(1800),
	68	fNEvents(0),
	69	fNAcceptedEvents(0)
	70	{
	71	//Default constructor. First we set the detector-type related constants.
	72	if (detectorType == kPhos) {
	73	fColumns = fgkPhosCols;
	74	fRows = fgkPhosRows;
	75	fModules = fgkPhosModules;
f4fc542c	76	fCaloString = "PHOS";
bd83f412	77	}
	78	else {
	79	//We'll just trust the enum to keep everything in line, so that if detectorType
	80	//isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
	81	//case, if someone intentionally gives another number
	82	fColumns = fgkEmCalCols;
	83	fRows = fgkEmCalRows;
	84	fModules = fgkEmCalModules;
f4fc542c	85	fCaloString = "EMCAL";
bd83f412	86	}
	87
	88	fDetType = detectorType;
	89
	90	// Set the number of points for each Amp vs. Time graph to 0
	91	memset(fNHighGain, 0, sizeof(fNHighGain));
	92	memset(fNLowGain, 0, sizeof(fNLowGain));
	93
	94	CreateGraphs(); // set up the TGraphs
	95
	96	// init TProfiles to NULL=0 also
	97	memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
	98	memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
	99	}
	100
	101	// dtor
	102	//_____________________________________________________________________
	103	AliCaloCalibSignal::~AliCaloCalibSignal()
	104	{
	105	ClearObjects();
	106	}
	107
	108	//_____________________________________________________________________
	109	void AliCaloCalibSignal::ClearObjects()
	110	{
	111	// delete what was created in the ctor (TGraphs), and possible later (TProfiles)
	112	for (int i=0; i<fgkMaxTowers; i++) {
	113	if ( fGraphAmpVsTimeHighGain[i] ) { delete fGraphAmpVsTimeHighGain[i]; }
	114	if ( fGraphAmpVsTimeLowGain[i] ) { delete fGraphAmpVsTimeLowGain[i]; }
	115	if ( fProfAmpVsTimeHighGain[i] ) { delete fProfAmpVsTimeHighGain[i]; }
	116	if ( fProfAmpVsTimeLowGain[i] ) { delete fProfAmpVsTimeLowGain[i]; }
	117	}
	118	// set pointers
	119	memset(fGraphAmpVsTimeHighGain, 0, sizeof(fGraphAmpVsTimeHighGain));
	120	memset(fGraphAmpVsTimeLowGain, 0, sizeof(fGraphAmpVsTimeLowGain));
	121	memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
	122	memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
	123
	124	return;
	125	}
	126
	127	// copy ctor
	128	//_____________________________________________________________________
	129	AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
	130	TObject(sig),
	131	fDetType(sig.GetDetectorType()),
	132	fColumns(sig.GetColumns()),
	133	fRows(sig.GetRows()),
	134	fModules(sig.GetModules()),
f4fc542c	135	fCaloString(sig.GetCaloString()),
f4fc542c	136	fMapping(NULL), //! note that we are not copying the map info
bd83f412	137	fRunNumber(sig.GetRunNumber()),
	138	fStartTime(sig.GetStartTime()),
	139	fAmpCut(sig.GetAmpCut()),
	140	fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
	141	fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
	142	fHour(sig.GetHour()),
	143	fLatestHour(sig.GetLatestHour()),
	144	fUseAverage(sig.GetUseAverage()),
	145	fSecInAverage(sig.GetSecInAverage()),
	146	fNEvents(sig.GetNEvents()),
	147	fNAcceptedEvents(sig.GetNAcceptedEvents())
	148	{
	149	// also the TGraph contents
	150	AddInfo(&sig);
	151	}
	152
	153	// assignment operator; use copy ctor to make life easy..
	154	//_____________________________________________________________________
	155	AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
	156	{
	157	// assignment operator; use copy ctor
	158	if (&source == this) return *this;
	159
	160	new (this) AliCaloCalibSignal(source);
	161	return *this;
	162	}
	163
	164	//_____________________________________________________________________
	165	void AliCaloCalibSignal::CreateGraphs()
	166	{
	167	//Then, loop for the requested number of modules
	168	TString title, name;
	169	for (int i = 0; i < fModules; i++) {
	170
	171	// Amplitude vs. Time graph for each channel
	172	for(int ic=0;ic < fColumns;ic++){
	173	for(int ir=0;ir < fRows;ir++){
	174
	175	int id = GetTowerNum(i, ic, ir);
	176
	177	// high gain graph
	178	name = "fGraphAmpVsTimeHighGain_"; name += i;
	179	name += "_"; name += ic;
	180	name += "_"; name += ir;
	181	title = "Amp vs. Time High Gain Mod "; title += i;
	182	title += " Col "; title += ic;
	183	title += " Row "; title += ir;
	184
	185	fGraphAmpVsTimeHighGain[id] = new TGraph();
	186	fGraphAmpVsTimeHighGain[id]->SetName(name);
	187	fGraphAmpVsTimeHighGain[id]->SetTitle(title);
	188	fGraphAmpVsTimeHighGain[id]->SetMarkerStyle(20);
	189
	190	// Low Gain
	191	name = "fGraphAmpVsTimeLowGain_"; name += i;
	192	name += "_"; name += ic;
	193	name += "_"; name += ir;
	194	title = "Amp vs. Time Low Gain Mod "; title += i;
	195	title += " Col "; title += ic;
	196	title += " Row "; title += ir;
	197
	198	fGraphAmpVsTimeLowGain[id] = new TGraph();
	199	fGraphAmpVsTimeLowGain[id]->SetName(name);
	200	fGraphAmpVsTimeLowGain[id]->SetTitle(title);
201	fGraphAmpVsTimeLowGain[id]->SetMarkerStyle(20);
202
203	}
204	}
205
206	}//end for nModules
207	}
208
209	//_____________________________________________________________________
210	void AliCaloCalibSignal::Reset()
211	{
212	Zero(); // set all counters to 0
213	ClearObjects(); // delete previous TGraphs and TProfiles
214	CreateGraphs(); // and create some new ones
215	return;
216	}
217
218	//_____________________________________________________________________
219	void AliCaloCalibSignal::Zero()
220	{
221	// set all counters to 0; not cuts etc.though
222	fHour = 0;
223	fLatestHour = 0;
224	fNEvents = 0;
225	fNAcceptedEvents = 0;
226	return;
227	}
228
229	//_____________________________________________________________________
230	Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
231	{
232	int nAbove = 0;
233
234	int TowerNum = 0;
235	for (int i = 0; i<fModules; i++) {
236	for (int j = 0; j<fColumns; j++) {
237	for (int k = 0; k<fRows; k++) {
238	TowerNum = GetTowerNum(i,j,k);
239	if (AmpVal[TowerNum] > fAmpCut) {
240	nAbove++;
241	}
242	}
243	}
244	}
245
246	double fraction = (1.0*nAbove) / nTotChan;
247
248	if (fraction > fReqFractionAboveAmpCutVal) {
249	return true;
250	}
251	else return false;
252	}
253
254	//_____________________________________________________________________
255	Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
256	{
257	// just do this for the basic graphs/profiles that get filled in ProcessEvent
258	// may not have data for all channels, but let's just Add everything..
4fb34e09	259	// Note: this method will run into problems with TProfile adding if the binning of
	260	// the local profiles is not the same as those provided by the argument *sig..
	261	int numGraphPoints = 0;
	262	int id = 0;
	263	int ip = 0;
bd83f412	264	for (int i = 0; i < fModules; i++) {
	265	for (int j = 0; j < fColumns; j++) {
	266	for (int k = 0; k < fRows; k++) {
	267
4fb34e09	268	id = GetTowerNum(i,j,k);
	269
	270	if(fUseAverage){ // add to Profiles
	271	if (sig->GetProfAmpVsTimeHighGain(id)) {
	272	GetProfAmpVsTimeHighGain(id)->Add(sig->GetProfAmpVsTimeHighGain(id));
	273	}
	274	if (sig->GetProfAmpVsTimeLowGain(id)) {
	275	GetProfAmpVsTimeLowGain(id)->Add(sig->GetProfAmpVsTimeLowGain(id));
	276	}
bd83f412	277	}
4fb34e09	278	else{ // add to Graphs
	279	// high gain
	280	numGraphPoints= sig->GetGraphAmpVsTimeHighGain(id)->GetN();
	281	if (numGraphPoints > 0) {
	282	// get the values
	283	double *graphX = sig->GetGraphAmpVsTimeHighGain(id)->GetX();
	284	double *graphY = sig->GetGraphAmpVsTimeHighGain(id)->GetY();
	285	for(ip=0; ip < numGraphPoints; ip++){
	286	fGraphAmpVsTimeHighGain[id]->SetPoint(fNHighGain[id]++,graphX[ip],graphY[ip]);
	287	}
	288	}
	289	// low gain
	290	numGraphPoints= sig->GetGraphAmpVsTimeLowGain(id)->GetN();
	291	if (numGraphPoints > 0) {
	292	// get the values
	293	double *graphX = sig->GetGraphAmpVsTimeLowGain(id)->GetX();
	294	double *graphY = sig->GetGraphAmpVsTimeLowGain(id)->GetY();
	295	for(ip=0; ip < numGraphPoints; ip++){
	296	fGraphAmpVsTimeLowGain[id]->SetPoint(fNLowGain[id]++,graphX[ip],graphY[ip]);
	297	}
	298	}
	299
bd83f412	300	}
	301
	302	}//end for nModules
	303	}//end for nColumns
	304	}//end for nRows
	305
	306	return kTRUE;//We succesfully added info from the supplied object
	307	}
	308
f4fc542c	309	//_____________________________________________________________________
	310	Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
	311	{
	312	// if fMapping is NULL the rawstream will crate its own mapping
	313	AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
	314
	315	return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
	316	}
bd83f412	317
	318	//_____________________________________________________________________
	319	Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream in, AliRawEventHeaderBase aliHeader)
	320	{
	321	// Method to process=analyze one event in the data stream
	322	if (!in) return kFALSE; //Return right away if there's a null pointer
	323
	324	fNEvents++; // one more event
	325
	326	// PHOS has more towers than EMCAL, so use PHOS numbers to set array sizes
	327	int AmpValHighGain[fgkMaxTowers];
	328	int AmpValLowGain[fgkMaxTowers];
	329
	330	memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
	331	memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
	332
8cdb1ddb	333	int sample, isample = 0; //The sample temp, and the sample number in current event.
bd83f412	334	int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal
	335	int gain = 0;
	336
	337	// Number of Low and High gain channels for this event:
	338	int nLowChan = 0;
	339	int nHighChan = 0;
	340
	341	int TowerNum = 0; // array index for TGraphs etc.
	342
	343	// loop first to get the fraction of channels with amplitudes above cut
	344	while (in->Next()) {
	345	sample = in->GetSignal(); //Get the adc signal
	346	if (sample < min) min = sample;
	347	if (sample > max) max = sample;
8cdb1ddb	348	isample++;
8cdb1ddb	349	if ( isample >= in->GetTimeLength()) {
bd83f412	350	//If we're here then we're done with this tower
	351	gain = 1 - in->IsLowGain();
	352
	353	int arrayPos = in->GetModule(); //The modules are numbered starting from 0
	354	if (arrayPos >= fModules) {
	355	//TODO: return an error message, if appopriate (perhaps if debug>0?)
	356	return kFALSE;
	357	}
	358
	359	//Debug
	360	if (arrayPos < 0 \|\| arrayPos >= fModules) {
	361	printf("Oh no: arrayPos = %i.\n", arrayPos);
	362	}
	363
	364	// get tower number for AmpVal array
	365	TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
	366
	367	if (gain == 0) {
	368	// fill amplitude into the array
	369	AmpValLowGain[TowerNum] = max - min;
	370	nLowChan++;
	371	}
	372	else if (gain==1) {//fill the high gain ones
	373	// fill amplitude into the array
	374	AmpValHighGain[TowerNum] = max - min;
	375	nHighChan++;
	376	}//end if gain
	377
	378
	379	max = fgkSampleMin; min = fgkSampleMax;
8cdb1ddb	380	isample = 0;
bd83f412	381
	382	}//End if end of tower
	383
	384	}//end while, of stream
	385
	386	// now check if it was a led event, only use high gain (that should be sufficient)
	387	if (fReqFractionAboveAmp) {
452729ca	388	bool ok = false;
	389	if (nHighChan > 0) {
	390	ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan);
	391	}
bd83f412	392	if (!ok) return false; // skip event
	393	}
	394
	395	fNAcceptedEvents++; // one more event accepted
	396
	397	if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
	398	fStartTime = aliHeader->Get("Timestamp");
	399	}
	400
	401	fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
	402	if (fLatestHour < fHour) {
	403	fLatestHour = fHour;
	404	}
	405
	406	// it is a led event, now fill graphs (maybe profiles later)
	407	for(int i=0;i<fModules;i++){
	408	for(int j=0;j<fColumns;j++){
	409	for(int k=0;k<fRows;k++){
	410
	411	TowerNum = GetTowerNum(i, j, k);
	412
	413	if(AmpValHighGain[TowerNum]) {
	414	fGraphAmpVsTimeHighGain[TowerNum]->SetPoint(fNHighGain[TowerNum]++,fHour,AmpValHighGain[TowerNum]);
	415	}
	416	if(AmpValLowGain[TowerNum]) {
	417	fGraphAmpVsTimeLowGain[TowerNum]->SetPoint(fNLowGain[TowerNum]++,fHour,AmpValLowGain[TowerNum]);
	418	}
	419	}
	420	}
	421	}
	422
	423	return kTRUE;
	424	}
	425
	426	//_____________________________________________________________________
	427	void AliCaloCalibSignal::CreateProfile(int imod, int ic, int ir, int towerId, int gain,
	428	int nbins, double min, double max)
	429	{ //! create/setup a TProfile
	430	TString title, name;
	431	if (gain == 0) {
	432	name = "fProfAmpVsTimeLowGain_";
	433	title = "Amp vs. Time Low Gain Mod ";
	434	}
	435	else if (gain == 1) {
	436	name = "fProfAmpVsTimeHighGain_";
	437	title = "Amp vs. Time High Gain Mod ";
	438	}
	439	name += imod;
	440	name += "_"; name += ic;
	441	name += "_"; name += ir;
	442	title += imod;
	443	title += " Col "; title += ic;
	444	title += " Row "; title += ir;
	445
	446	// use "s" option for RMS
	447	if (gain==0) {
	448	fProfAmpVsTimeLowGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
	449	}
	450	else if (gain==1) {
	451	fProfAmpVsTimeHighGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
	452	}
	453
	454	return;
	455	}
456	//_____________________________________________________________________
457	Bool_t AliCaloCalibSignal::Save(TString fileName, Bool_t saveEmptyGraphs)
458	{
459	//Saves all the histograms (or profiles, to be accurate) to the designated file
460
461	TFile destFile(fileName, "recreate");
462
463	if (destFile.IsZombie()) {
464	return kFALSE;
465	}
466
467	destFile.cd();
468
469	// setup variables for the TProfile plot
470	int numProfBins = 0;
471	double timeMin = 0;
472	double timeMax = 0;
473	if (fUseAverage) {
474	if (fSecInAverage > 0) {
475	numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
476	}
477	numProfBins += 2; // add extra buffer : first and last
478	double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
479	timeMin = - binSize;
480	timeMax = timeMin + numProfBins*binSize;
481	}
482
483	int numGraphPoints= 0;
484	int TowerNum = 0;
485	for (int i = 0; i < fModules; i++) {
486
487	for(int ic=0;ic < fColumns;ic++){
488	for(int ir=0;ir < fRows;ir++){
489
490	TowerNum = GetTowerNum(i, ic, ir);
491
492	// 1st: high gain
493	numGraphPoints= fGraphAmpVsTimeHighGain[TowerNum]->GetN();
494	if( numGraphPoints>0 \|\| saveEmptyGraphs) {
495
496	// average the graphs points over time if requested and put them in a profile plot
497	if(fUseAverage && numGraphPoints>0) {
498
499	// get the values
500	double *graphX = fGraphAmpVsTimeHighGain[TowerNum]->GetX();
501	double *graphY = fGraphAmpVsTimeHighGain[TowerNum]->GetY();
502
503	// create the TProfile: 1 is for High gain
504	CreateProfile(i, ic, ir, TowerNum, 1,
505	numProfBins, timeMin, timeMax);
506
507	// loop over graph points and fill profile
508	for(int ip=0; ip < numGraphPoints; ip++){
509	fProfAmpVsTimeHighGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
510	}
511
512	fProfAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
513	fProfAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
514	fProfAmpVsTimeHighGain[TowerNum]->Write();
515
516	}
517	else{
518	//otherwise, just save the graphs and forget the profiling
519	fGraphAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
520	fGraphAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
521	fGraphAmpVsTimeHighGain[TowerNum]->Write();
522	}
523
524	} // low gain graph info should be saved in one form or another
525
526	// 2nd: now go to the low gain case
527	numGraphPoints= fGraphAmpVsTimeLowGain[TowerNum]->GetN();
528	if( numGraphPoints>0 \|\| saveEmptyGraphs) {
529
530	// average the graphs points over time if requested and put them in a profile plot
531	if(fUseAverage && numGraphPoints>0) {
532
533	double *graphX = fGraphAmpVsTimeLowGain[TowerNum]->GetX();
534	double *graphY = fGraphAmpVsTimeLowGain[TowerNum]->GetY();
535
536	// create the TProfile: 0 is for Low gain
537	CreateProfile(i, ic, ir, TowerNum, 0,
538	numProfBins, timeMin, timeMax);
539
540	// loop over graph points and fill profile
541	for(int ip=0; ip < numGraphPoints; ip++){
542	fProfAmpVsTimeLowGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
543	}
544
545	fProfAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
546	fProfAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
547	fProfAmpVsTimeLowGain[TowerNum]->Write();
548
549	}
550	else{
551	//otherwise, just save the graphs and forget the profiling
552	fGraphAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
553	fGraphAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
554	fGraphAmpVsTimeLowGain[TowerNum]->Write();
555	}
556
557	} // low gain graph info should be saved in one form or another
558
559	} // fRows
560	} // fColumns
561
562	} // fModules
563	destFile.Close();
564
565	return kTRUE;
566	}