[u/mrichter/AliRoot.git] / FMD / AliFMDGainDA.cxx

/**************************************************************************
 * Copyright(c) 2004, 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.                  *
 **************************************************************************/

/** @file    AliFMDGainDA.cxx
    @author  Hans Hjersing Dalsgaard <canute@nbi.dk>
    @date    Mon Mar 13 13:46:05 2008
    @brief   Derived class for the pulse gain detector algorithm.
*/
//____________________________________________________________________
//
// This class contains the implementation of the gain detector
// algorithms (DA) for the FMD.  The data is collected in histograms
// that are reset for each pulse length after the mean and standard
// deviation are put into a TGraphErrors object. After a certain
// number of pulses (usually 8) the graph is fitted to a straight
// line. The gain is then slope of this line as it combines the known
// pulse and the response of the detector.
//
#include "AliFMDGainDA.h"
#include "AliFMDAltroMapping.h"
#include "AliFMDParameters.h"
#include "AliFMDCalibGain.h"
#include "AliFMDDigit.h"
#include "AliLog.h"
#include <TFile.h>
#include <TF1.h>
#include <TH1S.h>
#include <TGraphErrors.h>
#include <TDatime.h>
#include <TH2.h>

//_____________________________________________________________________
ClassImp(AliFMDGainDA)
#if 0 // Do not delete - here to let Emacs indent properly
;
#endif

//_____________________________________________________________________
AliFMDGainDA::AliFMDGainDA() 
  : AliFMDBaseDA(),
    fHighPulse(256), 
    fEventsPerChannel(10),
    fCurrentPulse(10),
    fCurrentChannel(10),
    fNumberOfStripsPerChip(128),
    fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
    fCurrentSummaryStrip(1),
    fGainFMD1i(0),
    fGainFMD2i(0),
    fGainFMD2o(0),
    fGainFMD3i(0),
    fGainFMD3o(0),
    fChi2FMD1i(0),
    fChi2FMD2i(0),
    fChi2FMD2o(0),
    fChi2FMD3i(0),
    fChi2FMD3o(0)
{
  // Constructor 
  // 
  // Parameters: 
  //   None
  fCurrentPulse.Reset(0);
  fCurrentChannel.Reset(0);
  fOutputFile.open("gains.csv");
  fDiagnosticsFilename = "diagnosticsGain.root";
}

//_____________________________________________________________________
AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) 
  : AliFMDBaseDA(gainDA),
    fHighPulse(gainDA.fHighPulse),
    fEventsPerChannel(gainDA.fEventsPerChannel),
    fCurrentPulse(gainDA.fCurrentPulse),
    fCurrentChannel(gainDA.fCurrentChannel),
    fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
    fSummaryGains(gainDA.fSummaryGains),
    fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip),
    fGainFMD1i(gainDA.fGainFMD1i),
    fGainFMD2i(gainDA.fGainFMD2i),
    fGainFMD2o(gainDA.fGainFMD2o),
    fGainFMD3i(gainDA.fGainFMD3i),
    fGainFMD3o(gainDA.fGainFMD3o),
    fChi2FMD1i(gainDA.fChi2FMD1i),
    fChi2FMD2i(gainDA.fChi2FMD2i),
    fChi2FMD2o(gainDA.fChi2FMD2o),
    fChi2FMD3i(gainDA.fChi2FMD3i),
    fChi2FMD3o(gainDA.fChi2FMD3o)
{  
  // Copy Constructor 
  // 
  // Parameters: 
  //   gainDA   Object to copy from 
  fCurrentPulse.Reset(0);
  fCurrentChannel.Reset(0);
}

//_____________________________________________________________________
AliFMDGainDA::~AliFMDGainDA() 
{
  // Destructor 
  // 
  // Parameters: 
  //   None
}

//_____________________________________________________________________
void AliFMDGainDA::Init() 
{
  // Initialize 
  // 
  // Parameters: 
  //   None
  Int_t nEventsRequired = 0;
  
  for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++) {
    Int_t nEvents = 0;
    if(fPulseSize.At(idx))
      nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
    fEventsPerChannel.AddAt(nEvents,idx);
    if(nEvents>nEventsRequired) 
      nEventsRequired = nEvents * fNumberOfStripsPerChip;    
  }
  SetRequiredEvents(nEventsRequired); 
}

//_____________________________________________________________________
void AliFMDGainDA::InitContainer(TDirectory* dir) 
{
  AliFMDBaseDA::InitContainer(dir);
  
  for(UShort_t det=1;det<=3;det++) {
    UShort_t sr = (det == 1 ? 1 : 0);
    for (UShort_t ir = sr; ir < 2; ir++) {
      Char_t   ring = (ir == 0 ? 'O' : 'I');
      UShort_t nsec = (ir == 0 ? 40  : 20);
      UShort_t nva  = (ir == 0 ? 2 : 4);
      for(UShort_t sec =0; sec < nsec;  sec++)  {
	TObjArray* sectorArray = GetSectorArray(det, ring, sec);
	TObjArray* cache = new TObjArray(nva);
	cache->SetName("Cache");
	cache->SetOwner();
	Int_t n = sectorArray->GetEntriesFast();
	sectorArray->AddAtAndExpand(cache, n);
	for(UShort_t va = 0; va < nva; va++) {
	  TH1S* hChannel = new TH1S(Form("FMD%d%c[%02d]_va%d",
					 det,ring,sec,va),
				    Form("FMD%d%c[%02d] VA%d Cache",
					 det,ring,sec,va),
				    1024,-.5,1023.5);
	  hChannel->SetDirectory(0);
	  cache->AddAtAndExpand(hChannel,va);
	}
      }
    }
  }
}

//_____________________________________________________________________
void AliFMDGainDA::AddChannelContainer(TObjArray* stripArray, 
				       UShort_t det  , 
				       Char_t   ring,  
				       UShort_t sec, 
				       UShort_t strip) 
{  
  // Make a channel container 
  // 
  // Parameters: 
  //  sectorArray    Sectors 
  //  det            Detector number
  //  ring           Ring identifier 
  //  sec            Sector number
  //  strip          Strip number

  AliFMDParameters* pars     = AliFMDParameters::Instance();
  UShort_t          board    = pars->GetAltroMap()->Sector2Board(ring, sec);
  Int_t             halfring = GetHalfringIndex(det,ring,board/16);
  Int_t             dPulse   = fPulseSize.At(halfring);
  Int_t             nPulses  = dPulse > 0 ? fHighPulse / dPulse : 0;

  TGraphErrors* gChannel  = new TGraphErrors(nPulses);
  gChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
  gChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC", 
			  det, ring, sec, strip));
  stripArray->AddAtAndExpand(gChannel,0);
}

//_____________________________________________________________________
void AliFMDGainDA::AddSectorSummary(TObjArray* sectorArray, 
				    UShort_t det, 
				    Char_t   ring, 
				    UShort_t sec, 
				    UShort_t nStr) 
{
  TH1F* summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]", 
					   det, ring, sec), 
			   nStr, -.5, nStr-.5);
  summary->SetXTitle("Strip");
  summary->SetYTitle("Gain [ADC/DAC]");
  summary->SetDirectory(0);

  Int_t n = sectorArray->GetEntriesFast();
  sectorArray->AddAtAndExpand(summary, n);
}

//_____________________________________________________________________
void AliFMDGainDA::FillChannels(AliFMDDigit* digit) 
{
  // Fill data into histogram
  // 
  // Parameters: 
  //   digit  Digit to get the data from

  UShort_t det   = digit->Detector();
  Char_t   ring  = digit->Ring();
  UShort_t sec   = digit->Sector();
  UShort_t strip = digit->Strip();
  
  //Strip is always seen as the first in a VA chip. All other strips are junk.
  //Strips are counted from zero on even sectors and from 511 on odd sectors...
   
  if((sec%2)     && ((strip+1) % fNumberOfStripsPerChip)) return;
  if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
  
  UShort_t vaChip   = strip / fNumberOfStripsPerChip; 
  TH1S*    hChannel = GetChannelHistogram(det, ring, sec, vaChip);
  hChannel->Fill(digit->Counts());
  UpdatePulseAndADC(det,ring,sec,strip);
}

//_____________________________________________________________________
void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
{
  // 
  // Create summary hists for FMD gains and chi2 of the fits
  //
  switch (det) { 
  case 1: 
    fGainFMD1i = MakeSummaryHistogram("gain", "Gains", det, ring);
    fChi2FMD1i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
    break;
  case 2:
    switch (ring) { 
    case 'I': case 'i':
      fGainFMD2i = MakeSummaryHistogram("gain", "Gains", det, ring);
      fChi2FMD2i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
      break;
    case 'O': case 'o':
      fGainFMD2o = MakeSummaryHistogram("gain", "Gains", det, ring);
      fChi2FMD2o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
      break;
    }
    break;
  case 3:
    switch (ring) { 
    case 'I': case 'i':
      fGainFMD3i = MakeSummaryHistogram("gain", "Gains", det, ring);
      fChi2FMD3i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
      break;
    case 'O': case 'o':
      fGainFMD3o = MakeSummaryHistogram("gain", "Gains", det, ring);
      fChi2FMD3o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
      break;
    }
    break;
  }
}

//_____________________________________________________________________
void AliFMDGainDA::Analyse(UShort_t det, 
			   Char_t   ring, 
			   UShort_t sec, 
			   UShort_t strip) 
{
  // Analyse result of a single strip
  // 
  // Parameters: 
  //  det            Detector number
  //  ring           Ring identifier 
  //  sec            Sector number
  //  strip          Strip number
  TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
  if(!grChannel->GetN()) {
    AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
                     det, ring , sec, strip));
    return;
  }
  TF1 fitFunc("fitFunc","pol1",-10,280); 
  fitFunc.SetParameters(100,3);
  grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
     
  Float_t gain    = -1;
  Float_t error   = -1; 
  Float_t chi2ndf = -1;
  if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
    gain    = fitFunc.GetParameter(1);
    error   = fitFunc.GetParError(1);
    if(fitFunc.GetNDF())
      chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
  }
  
  fOutputFile << det                         << ','
	      << ring                        << ','
	      << sec                         << ','
	      << strip                       << ','
	      << gain                        << ','
	      << error                       << ','
	      << chi2ndf                     <<"\n";
  
  //due to RCU trouble, first strips on VAs are excluded
  // if(strip%128 != 0) {
    
  fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
  fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
  
  fCurrentSummaryStrip++;

  TH2* hGain = 0;
  TH2* hChi2 = 0;
  switch (det) { 
  case 1: hGain = fGainFMD1i; hChi2 = fChi2FMD1i; break;
  case 2: 
    switch (ring) { 
    case 'I':  hGain = fGainFMD2i; hChi2 = fChi2FMD2i; break;
    case 'O':  hGain = fGainFMD2o; hChi2 = fChi2FMD2o; break;
    }
    break;
  case 3:
    switch (ring) { 
    case 'I':  hGain = fGainFMD3i; hChi2 = fChi2FMD3i; break;
    case 'O':  hGain = fGainFMD3o; hChi2 = fChi2FMD3o; break;
    }
    break;
  }
  if (hGain && hChi2) {
    Int_t bin = hGain->FindBin(sec, strip);
    hGain->SetBinContent(bin, gain);
    hGain->SetBinError(bin, error);
    hChi2->SetBinContent(bin, chi2ndf);
  }

    // }
  if(fSaveHistograms) {
    TH1F* summary = GetSectorSummary(det, ring, sec);
    summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
    summary->SetBinError(strip+1, fitFunc.GetParError(1));
  }  
}

//_____________________________________________________________________
void AliFMDGainDA::Terminate(TFile* diagFile)
{
  // End of file 
  // 
  // Parameters: 
  //   None
  if(diagFile) {
    diagFile->cd();
    fSummaryGains.Write();
  }
}

//_____________________________________________________________________
void AliFMDGainDA::WriteHeaderToFile() 
{
  // Write header to the output file
  // 
  // Parameters: 
  //   None
  AliFMDParameters* pars       = AliFMDParameters::Instance();
  fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
  TDatime now;
  fOutputFile << "# This file created from run # " << fRunno 
	      << " @ " << now.AsString() << std::endl;
  fOutputFile.write("# Detector, "
		    "Ring, "
		    "Sector, "
		    "Strip, "
		    "Gain, "
		    "Error, "
		    "Chi2/NDF \n",56);  
}

//_____________________________________________________________________
TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det, 
					Char_t   ring, 
					UShort_t sec, 
					UShort_t va) 
{
  // Get the current histogram of a single strip
  // 
  // Parameters: 
  //  det            Detector number
  //  ring           Ring identifier 
  //  sec            Sector number
  //  va             VA number
  TObjArray* secArray  = GetSectorArray(det, ring, sec);
  Int_t      n         = secArray->GetEntriesFast();
  TObjArray* cache     = static_cast<TObjArray*>(secArray->At(n-1));
  TH1S* hChannel       = static_cast<TH1S*>(cache->At(va));
  
  return hChannel;
}

//_____________________________________________________________________
TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det, 
				       Char_t   ring, 
				       UShort_t sec, 
				       UShort_t strip) 
{  
  // Get the graph of a single strip
  // 
  // Parameters: 
  //  det            Detector number
  //  ring           Ring identifier 
  //  sec            Sector number
  //  strip          Strip number
  TObjArray*    stripArray = GetStripArray(det, ring, sec, strip);
  TGraphErrors* hChannel   = static_cast<TGraphErrors*>(stripArray->At(0));
  
  return hChannel;
}

//_____________________________________________________________________
TH1F* AliFMDGainDA::GetSectorSummary(UShort_t det, 
				     Char_t   ring, 
				     UShort_t sec) 
{
  TObjArray* secArray    = GetSectorArray(det, ring, sec);
  Int_t      n           = secArray->GetEntriesFast();
  return static_cast<TH1F*>(secArray->At(n-2)); // Cache added later
}

//_____________________________________________________________________
void AliFMDGainDA::UpdatePulseAndADC(UShort_t det, 
				     Char_t ring, 
				     UShort_t sec, 
				     UShort_t strip) 
{
  // Go to next pulse size
  // 
  // Parameters: 
  //  det            Detector number
  //  ring           Ring identifier 
  //  sec            Sector number
  //  strip          Strip number
  
  AliFMDParameters* pars     = AliFMDParameters::Instance();
  UShort_t          board    = pars->GetAltroMap()->Sector2Board(ring, sec);
  Int_t             halfring = GetHalfringIndex(det,ring,board/16);
  
  if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
    return;
  
  if ((sec       % 2) && ((strip+1) % fNumberOfStripsPerChip)) return;
  if (((sec + 1) % 2) && (strip     % fNumberOfStripsPerChip)) return;
  
  if(((GetCurrentEvent()) % fPulseLength.At(halfring)) 
     && GetCurrentEvent() > 0) return;
     
  Int_t vaChip   = strip/fNumberOfStripsPerChip; 
  TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
  
  if(!hChannel->GetEntries()) {
    AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
		    det, ring , sec, strip));
    return;
  }
  Double_t mean      = hChannel->GetMean();
  Double_t rms       = hChannel->GetRMS();
  Double_t pulse     = (Double_t(fCurrentPulse.At(halfring)) 
			* fPulseSize.At(halfring));
  Int_t    firstBin  = hChannel->GetXaxis()->GetFirst();
  Int_t    lastBin   = hChannel->GetXaxis()->GetLast();
  hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
  
  mean               = hChannel->GetMean();
  rms                = hChannel->GetRMS();
  
  hChannel->GetXaxis()->SetRange(firstBin,lastBin);
  
  Int_t         channelOff = ((GetCurrentEvent()-1)
			      / ((fPulseLength.At(halfring)*fHighPulse)
				 / fPulseSize.At(halfring)));
  Int_t         channelNo  = (strip + ((sec % 2 == 1) ? -1 : 1) * channelOff);
  TGraphErrors* channel    = GetChannel(det,ring,sec,channelNo);

  channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
  channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
  
  if(fSaveHistograms) {
    TH1S* out = 
      static_cast<TH1S*>(hChannel->Clone(Form("FMD%d%c[%02d,%03d]_0x%02x",
					      det, ring, sec, channelNo, 
					      int(pulse))));
    out->SetTitle(Form("FMD%d%c[%02d,%03d] DAC=0x%02x (%3d)",
		       det, ring, sec, channelNo, int(pulse), int(pulse)));
    out->SetDirectory(0);
    TObjArray* arr = GetStripArray(det, ring, sec, channelNo);
    arr->AddAtAndExpand(out, fCurrentPulse.At(halfring)+1);
  }
    
  hChannel->Reset();
  
}

//_____________________________________________________________________
void AliFMDGainDA::ResetPulseAndUpdateChannel() 
{  
  // Reset all
  // 
  // Parameters: 
  //   None
  fCurrentPulse.Reset(0); 
}

//_____________________________________________________________________
void AliFMDGainDA::FinishEvent() 
{
  // End of event
  // 
  // Parameters: 
  //   None
  for(UShort_t det=1; det<=3;det++) {
    UShort_t firstring = (det == 1 ? 1 : 0);
    for(UShort_t iring = firstring; iring <=1;iring++) {
      Char_t ring = (iring == 1 ? 'I' : 'O');
      for(UShort_t board =0 ; board <=1; board++) {
	Int_t idx = GetHalfringIndex(det,ring,board);
	
	if(!fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
	  continue;
	if(GetCurrentEvent()>0 && 
	   ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
	  fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
	
	if(GetCurrentEvent()>0 && 
	   ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
	  fCurrentPulse.AddAt(0,idx);
      }
    }
  }
}
//_____________________________________________________________________
//
//EOF
//
Commit	Line	Data
f631f26d	1	/**************************************************************************
	2	* Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/** @file AliFMDGainDA.cxx
	17	@author Hans Hjersing Dalsgaard <canute@nbi.dk>
	18	@date Mon Mar 13 13:46:05 2008
	19	@brief Derived class for the pulse gain detector algorithm.
	20	*/
cf291b42	21	//____________________________________________________________________
	22	//
	23	// This class contains the implementation of the gain detector
	24	// algorithms (DA) for the FMD. The data is collected in histograms
	25	// that are reset for each pulse length after the mean and standard
	26	// deviation are put into a TGraphErrors object. After a certain
	27	// number of pulses (usually 8) the graph is fitted to a straight
	28	// line. The gain is then slope of this line as it combines the known
	29	// pulse and the response of the detector.
	30	//
f631f26d	31	#include "AliFMDGainDA.h"
b995fc28	32	#include "AliFMDAltroMapping.h"
09b6c804	33	#include "AliFMDParameters.h"
	34	#include "AliFMDCalibGain.h"
	35	#include "AliFMDDigit.h"
	36	#include "AliLog.h"
	37	#include <TFile.h>
	38	#include <TF1.h>
	39	#include <TH1S.h>
	40	#include <TGraphErrors.h>
6cf6e7a0	41	#include <TDatime.h>
2a082c96	42	#include <TH2.h>
f631f26d	43
	44	//_____________________________________________________________________
	45	ClassImp(AliFMDGainDA)
cf291b42	46	#if 0 // Do not delete - here to let Emacs indent properly
	47	;
	48	#endif
f631f26d	49
f631f26d	50	//_____________________________________________________________________
cf291b42	51	AliFMDGainDA::AliFMDGainDA()
cf291b42	52	: AliFMDBaseDA(),
cf291b42	53	fHighPulse(256),
22017a7e	54	fEventsPerChannel(10),
	55	fCurrentPulse(10),
	56	fCurrentChannel(10),
3bae5d02	57	fNumberOfStripsPerChip(128),
3bae5d02	58	fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
2a082c96	59	fCurrentSummaryStrip(1),
	60	fGainFMD1i(0),
	61	fGainFMD2i(0),
	62	fGainFMD2o(0),
	63	fGainFMD3i(0),
	64	fGainFMD3o(0),
	65	fChi2FMD1i(0),
	66	fChi2FMD2i(0),
	67	fChi2FMD2o(0),
	68	fChi2FMD3i(0),
	69	fChi2FMD3o(0)
f631f26d	70	{
6cf6e7a0	71	// Constructor
	72	//
	73	// Parameters:
	74	// None
427e8f99	75	fCurrentPulse.Reset(0);
427e8f99	76	fCurrentChannel.Reset(0);
f631f26d	77	fOutputFile.open("gains.csv");
a31ea3ce	78	fDiagnosticsFilename = "diagnosticsGain.root";
f631f26d	79	}
	80
	81	//_____________________________________________________________________
cf291b42	82	AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
2a082c96	83	: AliFMDBaseDA(gainDA),
2a082c96	84	fHighPulse(gainDA.fHighPulse),
	85	fEventsPerChannel(gainDA.fEventsPerChannel),
	86	fCurrentPulse(gainDA.fCurrentPulse),
	87	fCurrentChannel(gainDA.fCurrentChannel),
	88	fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
	89	fSummaryGains(gainDA.fSummaryGains),
	90	fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip),
	91	fGainFMD1i(gainDA.fGainFMD1i),
	92	fGainFMD2i(gainDA.fGainFMD2i),
	93	fGainFMD2o(gainDA.fGainFMD2o),
	94	fGainFMD3i(gainDA.fGainFMD3i),
	95	fGainFMD3o(gainDA.fGainFMD3o),
	96	fChi2FMD1i(gainDA.fChi2FMD1i),
	97	fChi2FMD2i(gainDA.fChi2FMD2i),
	98	fChi2FMD2o(gainDA.fChi2FMD2o),
	99	fChi2FMD3i(gainDA.fChi2FMD3i),
	100	fChi2FMD3o(gainDA.fChi2FMD3o)
cf291b42	101	{
6cf6e7a0	102	// Copy Constructor
	103	//
	104	// Parameters:
	105	// gainDA Object to copy from
427e8f99	106	fCurrentPulse.Reset(0);
427e8f99	107	fCurrentChannel.Reset(0);
f631f26d	108	}
	109
	110	//_____________________________________________________________________
cf291b42	111	AliFMDGainDA::~AliFMDGainDA()
cf291b42	112	{
6cf6e7a0	113	// Destructor
	114	//
	115	// Parameters:
	116	// None
f631f26d	117	}
	118
	119	//_____________________________________________________________________
cf291b42	120	void AliFMDGainDA::Init()
cf291b42	121	{
6cf6e7a0	122	// Initialize
	123	//
	124	// Parameters:
	125	// None
427e8f99	126	Int_t nEventsRequired = 0;
22017a7e	127
6cf6e7a0	128	for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++) {
	129	Int_t nEvents = 0;
	130	if(fPulseSize.At(idx))
	131	nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
	132	fEventsPerChannel.AddAt(nEvents,idx);
	133	if(nEvents>nEventsRequired)
a31ea3ce	134	nEventsRequired = nEvents * fNumberOfStripsPerChip;
6cf6e7a0	135	}
427e8f99	136	SetRequiredEvents(nEventsRequired);
a31ea3ce	137	}
	138
	139	//_____________________________________________________________________
	140	void AliFMDGainDA::InitContainer(TDirectory* dir)
	141	{
	142	AliFMDBaseDA::InitContainer(dir);
f631f26d	143
f631f26d	144	for(UShort_t det=1;det<=3;det++) {
a31ea3ce	145	UShort_t sr = (det == 1 ? 1 : 0);
a31ea3ce	146	for (UShort_t ir = sr; ir < 2; ir++) {
f631f26d	147	Char_t ring = (ir == 0 ? 'O' : 'I');
f631f26d	148	UShort_t nsec = (ir == 0 ? 40 : 20);
a31ea3ce	149	UShort_t nva = (ir == 0 ? 2 : 4);
f631f26d	150	for(UShort_t sec =0; sec < nsec; sec++) {
a31ea3ce	151	TObjArray* sectorArray = GetSectorArray(det, ring, sec);
	152	TObjArray* cache = new TObjArray(nva);
	153	cache->SetName("Cache");
	154	cache->SetOwner();
	155	Int_t n = sectorArray->GetEntriesFast();
	156	sectorArray->AddAtAndExpand(cache, n);
	157	for(UShort_t va = 0; va < nva; va++) {
	158	TH1S* hChannel = new TH1S(Form("FMD%d%c[%02d]_va%d",
	159	det,ring,sec,va),
	160	Form("FMD%d%c[%02d] VA%d Cache",
	161	det,ring,sec,va),
	162	1024,-.5,1023.5);
f631f26d	163	hChannel->SetDirectory(0);
a31ea3ce	164	cache->AddAtAndExpand(hChannel,va);
f631f26d	165	}
	166	}
	167	}
	168	}
	169	}
	170
	171	//_____________________________________________________________________
a31ea3ce	172	void AliFMDGainDA::AddChannelContainer(TObjArray* stripArray,
cf291b42	173	UShort_t det ,
cf291b42	174	Char_t ring,
f631f26d	175	UShort_t sec,
cf291b42	176	UShort_t strip)
cf291b42	177	{
6cf6e7a0	178	// Make a channel container
	179	//
	180	// Parameters:
	181	// sectorArray Sectors
	182	// det Detector number
	183	// ring Ring identifier
	184	// sec Sector number
	185	// strip Strip number
a31ea3ce	186
	187	AliFMDParameters* pars = AliFMDParameters::Instance();
	188	UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
	189	Int_t halfring = GetHalfringIndex(det,ring,board/16);
	190	Int_t dPulse = fPulseSize.At(halfring);
	191	Int_t nPulses = dPulse > 0 ? fHighPulse / dPulse : 0;
	192
	193	TGraphErrors* gChannel = new TGraphErrors(nPulses);
	194	gChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
	195	gChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC",
cf291b42	196	det, ring, sec, strip));
a31ea3ce	197	stripArray->AddAtAndExpand(gChannel,0);
	198	}
	199
	200	//_____________________________________________________________________
	201	void AliFMDGainDA::AddSectorSummary(TObjArray* sectorArray,
	202	UShort_t det,
	203	Char_t ring,
	204	UShort_t sec,
	205	UShort_t nStr)
	206	{
	207	TH1F* summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]",
	208	det, ring, sec),
	209	nStr, -.5, nStr-.5);
	210	summary->SetXTitle("Strip");
	211	summary->SetYTitle("Gain [ADC/DAC]");
	212	summary->SetDirectory(0);
	213
	214	Int_t n = sectorArray->GetEntriesFast();
	215	sectorArray->AddAtAndExpand(summary, n);
f631f26d	216	}
	217
	218	//_____________________________________________________________________
6cf6e7a0	219	void AliFMDGainDA::FillChannels(AliFMDDigit* digit)
	220	{
	221	// Fill data into histogram
	222	//
	223	// Parameters:
	224	// digit Digit to get the data from
f631f26d	225
	226	UShort_t det = digit->Detector();
	227	Char_t ring = digit->Ring();
	228	UShort_t sec = digit->Sector();
	229	UShort_t strip = digit->Strip();
	230
95779ff5	231	//Strip is always seen as the first in a VA chip. All other strips are junk.
	232	//Strips are counted from zero on even sectors and from 511 on odd sectors...
	233
	234	if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
	235	if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
	236
a31ea3ce	237	UShort_t vaChip = strip / fNumberOfStripsPerChip;
a31ea3ce	238	TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
f631f26d	239	hChannel->Fill(digit->Counts());
	240	UpdatePulseAndADC(det,ring,sec,strip);
	241	}
	242
2a082c96	243	//_____________________________________________________________________
	244	void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
	245	{
8cc1cb63	246	//
	247	// Create summary hists for FMD gains and chi2 of the fits
	248	//
2a082c96	249	switch (det) {
	250	case 1:
	251	fGainFMD1i = MakeSummaryHistogram("gain", "Gains", det, ring);
	252	fChi2FMD1i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	253	break;
	254	case 2:
	255	switch (ring) {
	256	case 'I': case 'i':
	257	fGainFMD2i = MakeSummaryHistogram("gain", "Gains", det, ring);
	258	fChi2FMD2i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	259	break;
	260	case 'O': case 'o':
	261	fGainFMD2o = MakeSummaryHistogram("gain", "Gains", det, ring);
	262	fChi2FMD2o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	263	break;
	264	}
	265	break;
	266	case 3:
	267	switch (ring) {
	268	case 'I': case 'i':
	269	fGainFMD3i = MakeSummaryHistogram("gain", "Gains", det, ring);
	270	fChi2FMD3i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	271	break;
	272	case 'O': case 'o':
	273	fGainFMD3o = MakeSummaryHistogram("gain", "Gains", det, ring);
	274	fChi2FMD3o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	275	break;
	276	}
	277	break;
	278	}
	279	}
	280
f631f26d	281	//_____________________________________________________________________
f631f26d	282	void AliFMDGainDA::Analyse(UShort_t det,
e9c06036	283	Char_t ring,
f631f26d	284	UShort_t sec,
6cf6e7a0	285	UShort_t strip)
	286	{
	287	// Analyse result of a single strip
	288	//
	289	// Parameters:
	290	// det Detector number
	291	// ring Ring identifier
	292	// sec Sector number
	293	// strip Strip number
f631f26d	294	TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
f631f26d	295	if(!grChannel->GetN()) {
95779ff5	296	AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
95779ff5	297	det, ring , sec, strip));
f631f26d	298	return;
f631f26d	299	}
9c958f2b	300	TF1 fitFunc("fitFunc","pol1",-10,280);
9c958f2b	301	fitFunc.SetParameters(100,3);
e9c06036	302	grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
f7f0b643	303
4fa150a7	304	Float_t gain = -1;
	305	Float_t error = -1;
	306	Float_t chi2ndf = -1;
	307	if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
	308	gain = fitFunc.GetParameter(1);
	309	error = fitFunc.GetParError(1);
	310	if(fitFunc.GetNDF())
	311	chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
	312	}
	313
f7f0b643	314	fOutputFile << det << ','
	315	<< ring << ','
	316	<< sec << ','
	317	<< strip << ','
4fa150a7	318	<< gain << ','
4fa150a7	319	<< error << ','
f631f26d	320	<< chi2ndf <<"\n";
f631f26d	321
22017a7e	322	//due to RCU trouble, first strips on VAs are excluded
276b1261	323	// if(strip%128 != 0) {
22017a7e	324
276b1261	325	fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
	326	fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
	327
	328	fCurrentSummaryStrip++;
2a082c96	329
	330	TH2* hGain = 0;
	331	TH2* hChi2 = 0;
	332	switch (det) {
	333	case 1: hGain = fGainFMD1i; hChi2 = fChi2FMD1i; break;
	334	case 2:
	335	switch (ring) {
	336	case 'I': hGain = fGainFMD2i; hChi2 = fChi2FMD2i; break;
	337	case 'O': hGain = fGainFMD2o; hChi2 = fChi2FMD2o; break;
	338	}
	339	break;
	340	case 3:
	341	switch (ring) {
	342	case 'I': hGain = fGainFMD3i; hChi2 = fChi2FMD3i; break;
	343	case 'O': hGain = fGainFMD3o; hChi2 = fChi2FMD3o; break;
	344	}
	345	break;
	346	}
	347	if (hGain && hChi2) {
	348	Int_t bin = hGain->FindBin(sec, strip);
	349	hGain->SetBinContent(bin, gain);
	350	hGain->SetBinError(bin, error);
	351	hChi2->SetBinContent(bin, chi2ndf);
	352	}
	353
276b1261	354	// }
f631f26d	355	if(fSaveHistograms) {
a31ea3ce	356	TH1F* summary = GetSectorSummary(det, ring, sec);
e9c06036	357	summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
e9c06036	358	summary->SetBinError(strip+1, fitFunc.GetParError(1));
e9c06036	359	}
f631f26d	360	}
	361
	362	//_____________________________________________________________________
3bae5d02	363	void AliFMDGainDA::Terminate(TFile* diagFile)
3bae5d02	364	{
6cf6e7a0	365	// End of file
	366	//
	367	// Parameters:
	368	// None
7bce699b	369	if(diagFile) {
	370	diagFile->cd();
	371	fSummaryGains.Write();
	372	}
3bae5d02	373	}
	374
	375	//_____________________________________________________________________
e9c06036	376	void AliFMDGainDA::WriteHeaderToFile()
e9c06036	377	{
6cf6e7a0	378	// Write header to the output file
	379	//
	380	// Parameters:
	381	// None
9c958f2b	382	AliFMDParameters* pars = AliFMDParameters::Instance();
9c958f2b	383	fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
6cf6e7a0	384	TDatime now;
	385	fOutputFile << "# This file created from run # " << fRunno
	386	<< " @ " << now.AsString() << std::endl;
f7f0b643	387	fOutputFile.write("# Detector, "
	388	"Ring, "
	389	"Sector, "
e9c06036	390	"Strip, "
	391	"Gain, "
	392	"Error, "
a31ea3ce	393	"Chi2/NDF \n",56);
f631f26d	394	}
	395
	396	//_____________________________________________________________________
	397	TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
e9c06036	398	Char_t ring,
f631f26d	399	UShort_t sec,
a31ea3ce	400	UShort_t va)
e9c06036	401	{
6cf6e7a0	402	// Get the current histogram of a single strip
	403	//
	404	// Parameters:
	405	// det Detector number
	406	// ring Ring identifier
	407	// sec Sector number
a31ea3ce	408	// va VA number
	409	TObjArray* secArray = GetSectorArray(det, ring, sec);
	410	Int_t n = secArray->GetEntriesFast();
	411	TObjArray* cache = static_cast<TObjArray*>(secArray->At(n-1));
	412	TH1S* hChannel = static_cast<TH1S*>(cache->At(va));
f631f26d	413
	414	return hChannel;
	415	}
	416
	417	//_____________________________________________________________________
cf291b42	418	TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
e9c06036	419	Char_t ring,
cf291b42	420	UShort_t sec,
e9c06036	421	UShort_t strip)
e9c06036	422	{
6cf6e7a0	423	// Get the graph of a single strip
	424	//
	425	// Parameters:
	426	// det Detector number
	427	// ring Ring identifier
	428	// sec Sector number
	429	// strip Strip number
a31ea3ce	430	TObjArray* stripArray = GetStripArray(det, ring, sec, strip);
a31ea3ce	431	TGraphErrors* hChannel = static_cast<TGraphErrors*>(stripArray->At(0));
f631f26d	432
	433	return hChannel;
	434	}
	435
a31ea3ce	436	//_____________________________________________________________________
	437	TH1F* AliFMDGainDA::GetSectorSummary(UShort_t det,
	438	Char_t ring,
	439	UShort_t sec)
	440	{
	441	TObjArray* secArray = GetSectorArray(det, ring, sec);
	442	Int_t n = secArray->GetEntriesFast();
	443	return static_cast<TH1F*>(secArray->At(n-2)); // Cache added later
	444	}
	445
f631f26d	446	//_____________________________________________________________________
cf291b42	447	void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
	448	Char_t ring,
	449	UShort_t sec,
	450	UShort_t strip)
	451	{
6cf6e7a0	452	// Go to next pulse size
	453	//
	454	// Parameters:
	455	// det Detector number
	456	// ring Ring identifier
	457	// sec Sector number
	458	// strip Strip number
f7f0b643	459
a31ea3ce	460	AliFMDParameters* pars = AliFMDParameters::Instance();
	461	UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
	462	Int_t halfring = GetHalfringIndex(det,ring,board/16);
22017a7e	463
f7f0b643	464	if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
aea2699c	465	return;
95779ff5	466
a31ea3ce	467	if ((sec % 2) && ((strip+1) % fNumberOfStripsPerChip)) return;
a31ea3ce	468	if (((sec + 1) % 2) && (strip % fNumberOfStripsPerChip)) return;
95779ff5	469
b995fc28	470	if(((GetCurrentEvent()) % fPulseLength.At(halfring))
b995fc28	471	&& GetCurrentEvent() > 0) return;
427e8f99	472
a31ea3ce	473	Int_t vaChip = strip/fNumberOfStripsPerChip;
cf291b42	474	TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
f631f26d	475
f631f26d	476	if(!hChannel->GetEntries()) {
cf291b42	477	AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
cf291b42	478	det, ring , sec, strip));
f631f26d	479	return;
	480	}
	481	Double_t mean = hChannel->GetMean();
	482	Double_t rms = hChannel->GetRMS();
b995fc28	483	Double_t pulse = (Double_t(fCurrentPulse.At(halfring))
b995fc28	484	* fPulseSize.At(halfring));
cf291b42	485	Int_t firstBin = hChannel->GetXaxis()->GetFirst();
cf291b42	486	Int_t lastBin = hChannel->GetXaxis()->GetLast();
9c958f2b	487	hChannel->GetXaxis()->SetRangeUser(mean-4rms,mean+4rms);
f631f26d	488
cf291b42	489	mean = hChannel->GetMean();
cf291b42	490	rms = hChannel->GetRMS();
efd39294	491
	492	hChannel->GetXaxis()->SetRange(firstBin,lastBin);
	493
a31ea3ce	494	Int_t channelOff = ((GetCurrentEvent()-1)
	495	/ ((fPulseLength.At(halfring)*fHighPulse)
	496	/ fPulseSize.At(halfring)));
	497	Int_t channelNo = (strip + ((sec % 2 == 1) ? -1 : 1) * channelOff);
	498	TGraphErrors* channel = GetChannel(det,ring,sec,channelNo);
	499
427e8f99	500	channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
427e8f99	501	channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
f631f26d	502
f631f26d	503	if(fSaveHistograms) {
a31ea3ce	504	TH1S* out =
	505	static_cast<TH1S*>(hChannel->Clone(Form("FMD%d%c[%02d,%03d]_0x%02x",
	506	det, ring, sec, channelNo,
	507	int(pulse))));
	508	out->SetTitle(Form("FMD%d%c[%02d,%03d] DAC=0x%02x (%3d)",
	509	det, ring, sec, channelNo, int(pulse), int(pulse)));
	510	out->SetDirectory(0);
	511	TObjArray* arr = GetStripArray(det, ring, sec, channelNo);
	512	arr->AddAtAndExpand(out, fCurrentPulse.At(halfring)+1);
f631f26d	513	}
22017a7e	514
f631f26d	515	hChannel->Reset();
	516
	517	}
	518
	519	//_____________________________________________________________________
cf291b42	520	void AliFMDGainDA::ResetPulseAndUpdateChannel()
cf291b42	521	{
6cf6e7a0	522	// Reset all
	523	//
	524	// Parameters:
	525	// None
427e8f99	526	fCurrentPulse.Reset(0);
f631f26d	527	}
	528
	529	//_____________________________________________________________________
cf291b42	530	void AliFMDGainDA::FinishEvent()
cf291b42	531	{
6cf6e7a0	532	// End of event
	533	//
	534	// Parameters:
	535	// None
22017a7e	536	for(UShort_t det=1; det<=3;det++) {
	537	UShort_t firstring = (det == 1 ? 1 : 0);
	538	for(UShort_t iring = firstring; iring <=1;iring++) {
	539	Char_t ring = (iring == 1 ? 'I' : 'O');
	540	for(UShort_t board =0 ; board <=1; board++) {
	541	Int_t idx = GetHalfringIndex(det,ring,board);
	542
6cf6e7a0	543	if(!fPulseLength.At(idx) \|\| !fEventsPerChannel.At(idx))
22017a7e	544	continue;
6cf6e7a0	545	if(GetCurrentEvent()>0 &&
6cf6e7a0	546	((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
22017a7e	547	fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
22017a7e	548
6cf6e7a0	549	if(GetCurrentEvent()>0 &&
6cf6e7a0	550	((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
22017a7e	551	fCurrentPulse.AddAt(0,idx);
	552	}
	553	}
427e8f99	554	}
f631f26d	555	}
	556	//_____________________________________________________________________
	557	//
	558	//EOF
	559	//