[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 <iostream>
// #include <fstream>
#include "AliLog.h"
#include "TF1.h"
#include "TH1S.h"
// #include "TMath.h"
#include "TGraphErrors.h"
#include "AliFMDParameters.h"
#include "AliFMDAltroMapping.h"
#include <TDatime.h>
#include <TH2.h>

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

//_____________________________________________________________________
AliFMDGainDA::AliFMDGainDA() 
  : AliFMDBaseDA(),
    fGainArray(),
    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");
  fGainArray.SetOwner(); 
}

//_____________________________________________________________________
AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) 
  : AliFMDBaseDA(gainDA),
    fGainArray(gainDA.fGainArray),
    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); 
  
  TObjArray* detArray;
  TObjArray* ringArray;
  TObjArray* sectorArray;
  
  for(UShort_t det=1;det<=3;det++) {
    detArray = new TObjArray();
    detArray->SetOwner();
    fGainArray.AddAtAndExpand(detArray,det);
    for (UShort_t ir = 0; ir < 2; ir++) {
      Char_t   ring = (ir == 0 ? 'O' : 'I');
      UShort_t nsec = (ir == 0 ? 40  : 20);
      UShort_t nstr = (ir == 0 ? 2 : 4);
      ringArray = new TObjArray();
      ringArray->SetOwner();
      detArray->AddAtAndExpand(ringArray,ir);
      for(UShort_t sec =0; sec < nsec;  sec++)  {
	sectorArray = new TObjArray();
	sectorArray->SetOwner();
	ringArray->AddAtAndExpand(sectorArray,sec);
	for(UShort_t strip = 0; strip < nstr; strip++) {
	  TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
				    Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
				    1024,0,1023);
	  hChannel->SetDirectory(0);
	  sectorArray->AddAtAndExpand(hChannel,strip);
	}
      }
    }
  }
}

//_____________________________________________________________________
void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray, 
				       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
  TGraphErrors* hChannel  = new TGraphErrors();
  hChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
  hChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC", 
			  det, ring, sec, strip));
  sectorArray->AddAtAndExpand(hChannel,strip);
}

//_____________________________________________________________________
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;
  
  Int_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)
{
  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) {
    gDirectory->cd(GetSectorPath(det,ring, sec, kTRUE));
    
    TH1F* summary = dynamic_cast<TH1F*>(gDirectory->Get("Summary"));
    if (!summary) { 
      Int_t nStr = (ring == 'I' ? 512 : 256);
      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(gDirectory);
    }
    summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
    summary->SetBinError(strip+1, fitFunc.GetParError(1));
    
    gDirectory->cd(GetStripPath(det,ring,sec,strip, kTRUE));
    grChannel->SetName(Form("FMD%d%c[%02d,%03d]",det,ring,sec,strip));
    // grChannel->SetDirectory(gDirectory);
    grChannel->Write();
    // grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
  }  
}

//_____________________________________________________________________
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 strip) 
{
  // Get the current histogram of a single strip
  // 
  // Parameters: 
  //  det            Detector number
  //  ring           Ring identifier 
  //  sec            Sector number
  //  strip          Strip number
  
  UShort_t  lRing = 1;
  if(ring == 'O')
    lRing = 0;
  
  
  TObjArray* detArray  = static_cast<TObjArray*>(fGainArray.At(det));
  TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(lRing));
  TObjArray* secArray  = static_cast<TObjArray*>(ringArray->At(sec));
  TH1S* hChannel       = static_cast<TH1S*>(secArray->At(strip));
  
  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
  UShort_t      iring     = (ring == 'O' ? 0 : 1);
  TObjArray*    detArray  = static_cast<TObjArray*>(fDetectorArray.At(det));
  TObjArray*    ringArray = static_cast<TObjArray*>(detArray->At(iring));
  TObjArray*    secArray  = static_cast<TObjArray*>(ringArray->At(sec));
  TGraphErrors* hChannel  = static_cast<TGraphErrors*>(secArray->At(strip));
  
  return hChannel;
}

//_____________________________________________________________________
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();
  // UInt_t ddl, board,chip,ch;
  UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
  // pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,ch);
  /// pars->GetAltroMap()->Strip2Channel(
  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    channelNumber      = (strip + 
				 (GetCurrentEvent()-1)
				 / ((fPulseLength.At(halfring)*fHighPulse)
				    / fPulseSize.At(halfring))); 
  if(sec%2)
    channelNumber      = (strip - 
			  (GetCurrentEvent()-1)
			  / ((fPulseLength.At(halfring)*fHighPulse)
			     / fPulseSize.At(halfring))); 
  
  TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
  
  channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
  channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
  
  if(fSaveHistograms) {
    gDirectory->cd(GetStripPath(det,ring,sec,channelNumber));
    hChannel->Write(Form("%s_pulse_%03d",hChannel->GetName(),(Int_t)pulse));
    
  }
    
  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"
6cf6e7a0	32	// #include <iostream>
6cf6e7a0	33	// #include <fstream>
f631f26d	34	#include "AliLog.h"
f631f26d	35	#include "TF1.h"
6cf6e7a0	36	#include "TH1S.h"
6cf6e7a0	37	// #include "TMath.h"
f631f26d	38	#include "TGraphErrors.h"
f631f26d	39	#include "AliFMDParameters.h"
b995fc28	40	#include "AliFMDAltroMapping.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()
	52	: AliFMDBaseDA(),
	53	fGainArray(),
cf291b42	54	fHighPulse(256),
22017a7e	55	fEventsPerChannel(10),
	56	fCurrentPulse(10),
	57	fCurrentChannel(10),
3bae5d02	58	fNumberOfStripsPerChip(128),
3bae5d02	59	fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
2a082c96	60	fCurrentSummaryStrip(1),
	61	fGainFMD1i(0),
	62	fGainFMD2i(0),
	63	fGainFMD2o(0),
	64	fGainFMD3i(0),
	65	fGainFMD3o(0),
	66	fChi2FMD1i(0),
	67	fChi2FMD2i(0),
	68	fChi2FMD2o(0),
	69	fChi2FMD3i(0),
	70	fChi2FMD3o(0)
f631f26d	71	{
6cf6e7a0	72	// Constructor
	73	//
	74	// Parameters:
	75	// None
427e8f99	76	fCurrentPulse.Reset(0);
427e8f99	77	fCurrentChannel.Reset(0);
f631f26d	78	fOutputFile.open("gains.csv");
cf291b42	79	fGainArray.SetOwner();
f631f26d	80	}
	81
	82	//_____________________________________________________________________
cf291b42	83	AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
2a082c96	84	: AliFMDBaseDA(gainDA),
	85	fGainArray(gainDA.fGainArray),
	86	fHighPulse(gainDA.fHighPulse),
	87	fEventsPerChannel(gainDA.fEventsPerChannel),
	88	fCurrentPulse(gainDA.fCurrentPulse),
	89	fCurrentChannel(gainDA.fCurrentChannel),
	90	fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
	91	fSummaryGains(gainDA.fSummaryGains),
	92	fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip),
	93	fGainFMD1i(gainDA.fGainFMD1i),
	94	fGainFMD2i(gainDA.fGainFMD2i),
	95	fGainFMD2o(gainDA.fGainFMD2o),
	96	fGainFMD3i(gainDA.fGainFMD3i),
	97	fGainFMD3o(gainDA.fGainFMD3o),
	98	fChi2FMD1i(gainDA.fChi2FMD1i),
	99	fChi2FMD2i(gainDA.fChi2FMD2i),
	100	fChi2FMD2o(gainDA.fChi2FMD2o),
	101	fChi2FMD3i(gainDA.fChi2FMD3i),
	102	fChi2FMD3o(gainDA.fChi2FMD3o)
cf291b42	103	{
6cf6e7a0	104	// Copy Constructor
	105	//
	106	// Parameters:
	107	// gainDA Object to copy from
427e8f99	108	fCurrentPulse.Reset(0);
427e8f99	109	fCurrentChannel.Reset(0);
f631f26d	110	}
	111
	112	//_____________________________________________________________________
cf291b42	113	AliFMDGainDA::~AliFMDGainDA()
cf291b42	114	{
6cf6e7a0	115	// Destructor
	116	//
	117	// Parameters:
	118	// None
f631f26d	119	}
	120
	121	//_____________________________________________________________________
cf291b42	122	void AliFMDGainDA::Init()
cf291b42	123	{
6cf6e7a0	124	// Initialize
	125	//
	126	// Parameters:
	127	// None
427e8f99	128	Int_t nEventsRequired = 0;
22017a7e	129
6cf6e7a0	130	for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++) {
	131	Int_t nEvents = 0;
	132	if(fPulseSize.At(idx))
	133	nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
	134	fEventsPerChannel.AddAt(nEvents,idx);
	135	if(nEvents>nEventsRequired)
	136	nEventsRequired = nEvents * fNumberOfStripsPerChip;
	137
	138	}
427e8f99	139
427e8f99	140	SetRequiredEvents(nEventsRequired);
22017a7e	141
f631f26d	142	TObjArray* detArray;
	143	TObjArray* ringArray;
	144	TObjArray* sectorArray;
	145
	146	for(UShort_t det=1;det<=3;det++) {
	147	detArray = new TObjArray();
	148	detArray->SetOwner();
	149	fGainArray.AddAtAndExpand(detArray,det);
	150	for (UShort_t ir = 0; ir < 2; ir++) {
	151	Char_t ring = (ir == 0 ? 'O' : 'I');
	152	UShort_t nsec = (ir == 0 ? 40 : 20);
	153	UShort_t nstr = (ir == 0 ? 2 : 4);
	154	ringArray = new TObjArray();
	155	ringArray->SetOwner();
	156	detArray->AddAtAndExpand(ringArray,ir);
	157	for(UShort_t sec =0; sec < nsec; sec++) {
	158	sectorArray = new TObjArray();
	159	sectorArray->SetOwner();
	160	ringArray->AddAtAndExpand(sectorArray,sec);
	161	for(UShort_t strip = 0; strip < nstr; strip++) {
cf291b42	162	TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
	163	Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
	164	1024,0,1023);
f631f26d	165	hChannel->SetDirectory(0);
	166	sectorArray->AddAtAndExpand(hChannel,strip);
	167	}
	168	}
	169	}
	170	}
	171	}
	172
	173	//_____________________________________________________________________
	174	void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray,
cf291b42	175	UShort_t det ,
cf291b42	176	Char_t ring,
f631f26d	177	UShort_t sec,
cf291b42	178	UShort_t strip)
cf291b42	179	{
6cf6e7a0	180	// Make a channel container
	181	//
	182	// Parameters:
	183	// sectorArray Sectors
	184	// det Detector number
	185	// ring Ring identifier
	186	// sec Sector number
	187	// strip Strip number
f631f26d	188	TGraphErrors* hChannel = new TGraphErrors();
cf291b42	189	hChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
	190	hChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC",
	191	det, ring, sec, strip));
f631f26d	192	sectorArray->AddAtAndExpand(hChannel,strip);
	193	}
	194
	195	//_____________________________________________________________________
6cf6e7a0	196	void AliFMDGainDA::FillChannels(AliFMDDigit* digit)
	197	{
	198	// Fill data into histogram
	199	//
	200	// Parameters:
	201	// digit Digit to get the data from
f631f26d	202
	203	UShort_t det = digit->Detector();
	204	Char_t ring = digit->Ring();
	205	UShort_t sec = digit->Sector();
	206	UShort_t strip = digit->Strip();
	207
95779ff5	208	//Strip is always seen as the first in a VA chip. All other strips are junk.
	209	//Strips are counted from zero on even sectors and from 511 on odd sectors...
	210
	211	if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
	212	if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
	213
cf291b42	214	Int_t vaChip = strip / fNumberOfStripsPerChip;
cf291b42	215	TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
f631f26d	216	hChannel->Fill(digit->Counts());
	217	UpdatePulseAndADC(det,ring,sec,strip);
	218	}
	219
2a082c96	220	//_____________________________________________________________________
	221	void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
	222	{
	223	switch (det) {
	224	case 1:
	225	fGainFMD1i = MakeSummaryHistogram("gain", "Gains", det, ring);
	226	fChi2FMD1i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	227	break;
	228	case 2:
	229	switch (ring) {
	230	case 'I': case 'i':
	231	fGainFMD2i = MakeSummaryHistogram("gain", "Gains", det, ring);
	232	fChi2FMD2i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	233	break;
	234	case 'O': case 'o':
	235	fGainFMD2o = MakeSummaryHistogram("gain", "Gains", det, ring);
	236	fChi2FMD2o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	237	break;
	238	}
	239	break;
	240	case 3:
	241	switch (ring) {
	242	case 'I': case 'i':
	243	fGainFMD3i = MakeSummaryHistogram("gain", "Gains", det, ring);
	244	fChi2FMD3i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	245	break;
	246	case 'O': case 'o':
	247	fGainFMD3o = MakeSummaryHistogram("gain", "Gains", det, ring);
	248	fChi2FMD3o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
	249	break;
	250	}
	251	break;
	252	}
	253	}
	254
f631f26d	255	//_____________________________________________________________________
f631f26d	256	void AliFMDGainDA::Analyse(UShort_t det,
e9c06036	257	Char_t ring,
f631f26d	258	UShort_t sec,
6cf6e7a0	259	UShort_t strip)
	260	{
	261	// Analyse result of a single strip
	262	//
	263	// Parameters:
	264	// det Detector number
	265	// ring Ring identifier
	266	// sec Sector number
	267	// strip Strip number
f631f26d	268	TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
f631f26d	269	if(!grChannel->GetN()) {
95779ff5	270	AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
95779ff5	271	det, ring , sec, strip));
f631f26d	272	return;
f631f26d	273	}
9c958f2b	274	TF1 fitFunc("fitFunc","pol1",-10,280);
9c958f2b	275	fitFunc.SetParameters(100,3);
e9c06036	276	grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
f7f0b643	277
4fa150a7	278	Float_t gain = -1;
	279	Float_t error = -1;
	280	Float_t chi2ndf = -1;
	281	if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
	282	gain = fitFunc.GetParameter(1);
	283	error = fitFunc.GetParError(1);
	284	if(fitFunc.GetNDF())
	285	chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
	286	}
	287
f7f0b643	288	fOutputFile << det << ','
	289	<< ring << ','
	290	<< sec << ','
	291	<< strip << ','
4fa150a7	292	<< gain << ','
4fa150a7	293	<< error << ','
f631f26d	294	<< chi2ndf <<"\n";
f631f26d	295
22017a7e	296	//due to RCU trouble, first strips on VAs are excluded
276b1261	297	// if(strip%128 != 0) {
22017a7e	298
276b1261	299	fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
	300	fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
	301
	302	fCurrentSummaryStrip++;
2a082c96	303
	304	TH2* hGain = 0;
	305	TH2* hChi2 = 0;
	306	switch (det) {
	307	case 1: hGain = fGainFMD1i; hChi2 = fChi2FMD1i; break;
	308	case 2:
	309	switch (ring) {
	310	case 'I': hGain = fGainFMD2i; hChi2 = fChi2FMD2i; break;
	311	case 'O': hGain = fGainFMD2o; hChi2 = fChi2FMD2o; break;
	312	}
	313	break;
	314	case 3:
	315	switch (ring) {
	316	case 'I': hGain = fGainFMD3i; hChi2 = fChi2FMD3i; break;
	317	case 'O': hGain = fGainFMD3o; hChi2 = fChi2FMD3o; break;
	318	}
	319	break;
	320	}
	321	if (hGain && hChi2) {
	322	Int_t bin = hGain->FindBin(sec, strip);
	323	hGain->SetBinContent(bin, gain);
	324	hGain->SetBinError(bin, error);
	325	hChi2->SetBinContent(bin, chi2ndf);
	326	}
	327
276b1261	328	// }
f631f26d	329	if(fSaveHistograms) {
e9c06036	330	gDirectory->cd(GetSectorPath(det,ring, sec, kTRUE));
	331
	332	TH1F* summary = dynamic_cast<TH1F*>(gDirectory->Get("Summary"));
	333	if (!summary) {
	334	Int_t nStr = (ring == 'I' ? 512 : 256);
	335	summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]",
	336	det, ring, sec),
	337	nStr, -.5, nStr-.5);
	338	summary->SetXTitle("Strip");
	339	summary->SetYTitle("Gain [ADC/DAC]");
	340	summary->SetDirectory(gDirectory);
	341	}
	342	summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
	343	summary->SetBinError(strip+1, fitFunc.GetParError(1));
	344
	345	gDirectory->cd(GetStripPath(det,ring,sec,strip, kTRUE));
	346	grChannel->SetName(Form("FMD%d%c[%02d,%03d]",det,ring,sec,strip));
	347	// grChannel->SetDirectory(gDirectory);
	348	grChannel->Write();
	349	// grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
	350	}
f631f26d	351	}
	352
	353	//_____________________________________________________________________
3bae5d02	354	void AliFMDGainDA::Terminate(TFile* diagFile)
3bae5d02	355	{
6cf6e7a0	356	// End of file
	357	//
	358	// Parameters:
	359	// None
7bce699b	360	if(diagFile) {
	361	diagFile->cd();
	362	fSummaryGains.Write();
	363	}
3bae5d02	364	}
	365
	366	//_____________________________________________________________________
e9c06036	367	void AliFMDGainDA::WriteHeaderToFile()
e9c06036	368	{
6cf6e7a0	369	// Write header to the output file
	370	//
	371	// Parameters:
	372	// None
9c958f2b	373	AliFMDParameters* pars = AliFMDParameters::Instance();
9c958f2b	374	fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
6cf6e7a0	375	TDatime now;
	376	fOutputFile << "# This file created from run # " << fRunno
	377	<< " @ " << now.AsString() << std::endl;
f7f0b643	378	fOutputFile.write("# Detector, "
	379	"Ring, "
	380	"Sector, "
e9c06036	381	"Strip, "
	382	"Gain, "
	383	"Error, "
f7f0b643	384	"Chi2/NDF \n",56);
f631f26d	385
	386	}
	387
	388	//_____________________________________________________________________
	389	TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
e9c06036	390	Char_t ring,
f631f26d	391	UShort_t sec,
e9c06036	392	UShort_t strip)
e9c06036	393	{
6cf6e7a0	394	// Get the current histogram of a single strip
	395	//
	396	// Parameters:
	397	// det Detector number
	398	// ring Ring identifier
	399	// sec Sector number
	400	// strip Strip number
f631f26d	401
6cf6e7a0	402	UShort_t lRing = 1;
f631f26d	403	if(ring == 'O')
6cf6e7a0	404	lRing = 0;
f631f26d	405
	406
	407	TObjArray* detArray = static_cast<TObjArray*>(fGainArray.At(det));
6cf6e7a0	408	TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(lRing));
f631f26d	409	TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
	410	TH1S* hChannel = static_cast<TH1S*>(secArray->At(strip));
	411
	412	return hChannel;
	413	}
	414
	415	//_____________________________________________________________________
cf291b42	416	TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
e9c06036	417	Char_t ring,
cf291b42	418	UShort_t sec,
e9c06036	419	UShort_t strip)
e9c06036	420	{
6cf6e7a0	421	// Get the graph of a single strip
	422	//
	423	// Parameters:
	424	// det Detector number
	425	// ring Ring identifier
	426	// sec Sector number
	427	// strip Strip number
cf291b42	428	UShort_t iring = (ring == 'O' ? 0 : 1);
	429	TObjArray* detArray = static_cast<TObjArray*>(fDetectorArray.At(det));
	430	TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(iring));
	431	TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
	432	TGraphErrors* hChannel = static_cast<TGraphErrors*>(secArray->At(strip));
f631f26d	433
	434	return hChannel;
	435	}
	436
	437	//_____________________________________________________________________
cf291b42	438	void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
	439	Char_t ring,
	440	UShort_t sec,
	441	UShort_t strip)
	442	{
6cf6e7a0	443	// Go to next pulse size
	444	//
	445	// Parameters:
	446	// det Detector number
	447	// ring Ring identifier
	448	// sec Sector number
	449	// strip Strip number
f7f0b643	450
427e8f99	451	AliFMDParameters* pars = AliFMDParameters::Instance();
b995fc28	452	// UInt_t ddl, board,chip,ch;
	453	UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
	454	// pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,ch);
	455	/// pars->GetAltroMap()->Strip2Channel(
22017a7e	456	Int_t halfring = GetHalfringIndex(det,ring,board/16);
22017a7e	457
f7f0b643	458	if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
aea2699c	459	return;
95779ff5	460
	461	if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
	462
	463	if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
	464
b995fc28	465	if(((GetCurrentEvent()) % fPulseLength.At(halfring))
b995fc28	466	&& GetCurrentEvent() > 0) return;
427e8f99	467
cf291b42	468	Int_t vaChip = strip/fNumberOfStripsPerChip;
cf291b42	469	TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
f631f26d	470
f631f26d	471	if(!hChannel->GetEntries()) {
cf291b42	472	AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
cf291b42	473	det, ring , sec, strip));
f631f26d	474	return;
	475	}
	476	Double_t mean = hChannel->GetMean();
	477	Double_t rms = hChannel->GetRMS();
b995fc28	478	Double_t pulse = (Double_t(fCurrentPulse.At(halfring))
b995fc28	479	* fPulseSize.At(halfring));
cf291b42	480	Int_t firstBin = hChannel->GetXaxis()->GetFirst();
cf291b42	481	Int_t lastBin = hChannel->GetXaxis()->GetLast();
9c958f2b	482	hChannel->GetXaxis()->SetRangeUser(mean-4rms,mean+4rms);
f631f26d	483
cf291b42	484	mean = hChannel->GetMean();
cf291b42	485	rms = hChannel->GetRMS();
efd39294	486
	487	hChannel->GetXaxis()->SetRange(firstBin,lastBin);
	488
95779ff5	489	Int_t channelNumber = (strip +
	490	(GetCurrentEvent()-1)
	491	/ ((fPulseLength.At(halfring)*fHighPulse)
	492	/ fPulseSize.At(halfring)));
	493	if(sec%2)
	494	channelNumber = (strip -
	495	(GetCurrentEvent()-1)
	496	/ ((fPulseLength.At(halfring)*fHighPulse)
	497	/ fPulseSize.At(halfring)));
f631f26d	498
	499	TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
	500
427e8f99	501	channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
427e8f99	502	channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
f631f26d	503
f631f26d	504	if(fSaveHistograms) {
f7f0b643	505	gDirectory->cd(GetStripPath(det,ring,sec,channelNumber));
	506	hChannel->Write(Form("%s_pulse_%03d",hChannel->GetName(),(Int_t)pulse));
	507
f631f26d	508	}
22017a7e	509
f631f26d	510	hChannel->Reset();
	511
	512	}
	513
	514	//_____________________________________________________________________
cf291b42	515	void AliFMDGainDA::ResetPulseAndUpdateChannel()
cf291b42	516	{
6cf6e7a0	517	// Reset all
	518	//
	519	// Parameters:
	520	// None
427e8f99	521	fCurrentPulse.Reset(0);
f631f26d	522	}
	523
	524	//_____________________________________________________________________
cf291b42	525	void AliFMDGainDA::FinishEvent()
cf291b42	526	{
6cf6e7a0	527	// End of event
	528	//
	529	// Parameters:
	530	// None
22017a7e	531	for(UShort_t det=1; det<=3;det++) {
	532	UShort_t firstring = (det == 1 ? 1 : 0);
	533	for(UShort_t iring = firstring; iring <=1;iring++) {
	534	Char_t ring = (iring == 1 ? 'I' : 'O');
	535	for(UShort_t board =0 ; board <=1; board++) {
	536	Int_t idx = GetHalfringIndex(det,ring,board);
	537
6cf6e7a0	538	if(!fPulseLength.At(idx) \|\| !fEventsPerChannel.At(idx))
22017a7e	539	continue;
6cf6e7a0	540	if(GetCurrentEvent()>0 &&
6cf6e7a0	541	((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
22017a7e	542	fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
22017a7e	543
6cf6e7a0	544	if(GetCurrentEvent()>0 &&
6cf6e7a0	545	((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
22017a7e	546	fCurrentPulse.AddAt(0,idx);
	547	}
	548	}
427e8f99	549	}
f631f26d	550	}
	551	//_____________________________________________________________________
	552	//
	553	//EOF
	554	//