[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 "TH1.h"
#include "TMath.h"
#include "TGraphErrors.h"
#include "AliFMDParameters.h"

//_____________________________________________________________________
ClassImp(AliFMDGainDA)

//_____________________________________________________________________
AliFMDGainDA::AliFMDGainDA() : AliFMDBaseDA(),
  fGainArray(),
  fPulseSize(32),
  fHighPulse(256), 
  fPulseLength(100),
  fEventsPerChannel(0),
  fCurrentPulse(0),
  fCurrentChannel(0),
  fNumberOfStripsPerChip(128)
{
  fOutputFile.open("gains.csv");
  fGainArray.SetOwner();
  
}

//_____________________________________________________________________
AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) :  
  AliFMDBaseDA(gainDA),
  fGainArray(gainDA.fGainArray),
  fPulseSize(gainDA.fPulseSize),
  fHighPulse(gainDA.fHighPulse),
  fPulseLength(gainDA.fPulseLength),
  fEventsPerChannel(gainDA.fEventsPerChannel),
  fCurrentPulse(gainDA.fCurrentPulse),
  fCurrentChannel(gainDA.fCurrentChannel),
  fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip)
{
  
}

//_____________________________________________________________________
AliFMDGainDA::~AliFMDGainDA() {
  
 
}

//_____________________________________________________________________
void AliFMDGainDA::Init() {
    
  fEventsPerChannel = (fPulseLength*fHighPulse) / fPulseSize ;
  
  SetRequiredEvents(fEventsPerChannel*fNumberOfStripsPerChip); //8 pulser values * 128 strips * 100 samples
  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) {
  
  TGraphErrors* hChannel  = new TGraphErrors();
  sectorArray->AddAtAndExpand(hChannel,strip);
}

//_____________________________________________________________________
void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {

  UShort_t det   = digit->Detector();
  Char_t   ring  = digit->Ring();
  UShort_t sec   = digit->Sector();
  UShort_t strip = digit->Strip();
  
  
  if(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::Analyse(UShort_t det, 
			   Char_t ring, 
			   UShort_t sec, 
			   UShort_t strip) {
  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","Q","Q",0,fHighPulse);
  AliFMDParameters* pars = AliFMDParameters::Instance();
  UInt_t ddl, board,chip,channel;
  pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,channel);
    
  Float_t chi2ndf = 0;
  if(fitFunc.GetNDF())
    chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
  ddl = ddl + kBaseDDL;
  
  Int_t relStrip = strip%fNumberOfStripsPerChip;
  
  fOutputFile << ddl                         << ','
	      << board                       << ','
	      << chip                        << ','
	      << channel                     << ','
	      << relStrip                    << ','
	      << fitFunc.GetParameter(1)     << ','
	      << fitFunc.GetParError(1)      << ','
	      << chi2ndf                     <<"\n";
  
  
  if(fSaveHistograms) {
    
    gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,strip));
    grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
  }
  
 
}

//_____________________________________________________________________
void AliFMDGainDA::WriteHeaderToFile() {
  AliFMDParameters* pars       = AliFMDParameters::Instance();
  fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
  fOutputFile.write("# Rcu, Board, Chip, Channel, Strip, Gain, Error, Chi2/NDF \n",59);
  
}

//_____________________________________________________________________
TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det, 
					Char_t ring, 
					UShort_t sec, 
					UShort_t strip) {
  
  UShort_t  Ring = 1;
  if(ring == 'O')
    Ring = 0;
  
  
  TObjArray* detArray  = static_cast<TObjArray*>(fGainArray.At(det));
  TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
  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) {
  
  UShort_t  Ring = 1;
  if(ring == 'O')
    Ring = 0;
  
  
  TObjArray* detArray          = static_cast<TObjArray*>(fDetectorArray.At(det));
  TObjArray* ringArray         = static_cast<TObjArray*>(detArray->At(Ring));
  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) {
  
  if(strip%fNumberOfStripsPerChip) return;
  if(((GetCurrentEvent())%fPulseLength) && 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*fPulseSize;
  
  hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
  
  mean      = hChannel->GetMean();
  rms       = hChannel->GetRMS();
    
  Int_t channelNumber = strip + (GetCurrentEvent()-1)/800; 
  
  TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
  
  channel->SetPoint(fCurrentPulse,pulse,mean);
  channel->SetPointError(fCurrentPulse,0,rms);
  
  if(fSaveHistograms) {
    gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,channelNumber));
    hChannel->GetXaxis()->SetRange(0,1023);
    hChannel->Write(Form("hFMD%d%c_%d_%d_pulse_%d",det,ring,sec,channelNumber,fCurrentPulse));
  }
  
  hChannel->Reset();
  
}

//_____________________________________________________________________
void AliFMDGainDA::ResetPulseAndUpdateChannel() {
  
  fCurrentPulse   = 0;
  
}

//_____________________________________________________________________
void AliFMDGainDA::FinishEvent() {
  
  if(GetCurrentEvent()>0 && (GetCurrentEvent()%fPulseLength==0))
    fCurrentPulse++;
  
  if(GetCurrentEvent()>0 && (GetCurrentEvent())%fEventsPerChannel==0)
    fCurrentPulse = 0;


}
//_____________________________________________________________________
//
//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	*/
	21	// This class contains the implementation of the gain detector algorithms (DA) for the FMD.
	22	// The data is collected in histograms that are reset for each pulse length after the
	23	// mean and standard deviation are put into a TGraphErrors object. After a certain number of pulses
	24	// (usually 8) the graph is fitted to a straight line. The gain is then slope of this line as
	25	// it combines the known pulse and the response of the detector.
	26
	27
	28	#include "AliFMDGainDA.h"
	29	#include "iostream"
	30	#include "fstream"
	31	#include "AliLog.h"
	32	#include "TF1.h"
	33	#include "TH1.h"
	34	#include "TMath.h"
	35	#include "TGraphErrors.h"
	36	#include "AliFMDParameters.h"
	37
	38	//_____________________________________________________________________
	39	ClassImp(AliFMDGainDA)
	40
	41	//_____________________________________________________________________
	42	AliFMDGainDA::AliFMDGainDA() : AliFMDBaseDA(),
	43	fGainArray(),
	44	fPulseSize(32),
	45	fHighPulse(256),
	46	fPulseLength(100),
	47	fEventsPerChannel(0),
	48	fCurrentPulse(0),
	49	fCurrentChannel(0),
	50	fNumberOfStripsPerChip(128)
	51	{
	52	fOutputFile.open("gains.csv");
	53	fGainArray.SetOwner();
	54
	55	}
	56
	57	//_____________________________________________________________________
	58	AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) :
	59	AliFMDBaseDA(gainDA),
	60	fGainArray(gainDA.fGainArray),
	61	fPulseSize(gainDA.fPulseSize),
	62	fHighPulse(gainDA.fHighPulse),
	63	fPulseLength(gainDA.fPulseLength),
	64	fEventsPerChannel(gainDA.fEventsPerChannel),
65	fCurrentPulse(gainDA.fCurrentPulse),
66	fCurrentChannel(gainDA.fCurrentChannel),
67	fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip)
68	{
69
70	}
71
72	//_____________________________________________________________________
73	AliFMDGainDA::~AliFMDGainDA() {
74
75
76
77	}
78
79	//_____________________________________________________________________
80	void AliFMDGainDA::Init() {
f631f26d	81
	82	fEventsPerChannel = (fPulseLength*fHighPulse) / fPulseSize ;
	83
	84	SetRequiredEvents(fEventsPerChannelfNumberOfStripsPerChip); //8 pulser values 128 strips * 100 samples
	85	TObjArray* detArray;
	86	TObjArray* ringArray;
	87	TObjArray* sectorArray;
	88
	89	for(UShort_t det=1;det<=3;det++) {
	90	detArray = new TObjArray();
	91	detArray->SetOwner();
	92	fGainArray.AddAtAndExpand(detArray,det);
	93	for (UShort_t ir = 0; ir < 2; ir++) {
	94	Char_t ring = (ir == 0 ? 'O' : 'I');
	95	UShort_t nsec = (ir == 0 ? 40 : 20);
	96	UShort_t nstr = (ir == 0 ? 2 : 4);
	97	ringArray = new TObjArray();
	98	ringArray->SetOwner();
	99	detArray->AddAtAndExpand(ringArray,ir);
	100	for(UShort_t sec =0; sec < nsec; sec++) {
	101	sectorArray = new TObjArray();
	102	sectorArray->SetOwner();
	103	ringArray->AddAtAndExpand(sectorArray,sec);
	104	for(UShort_t strip = 0; strip < nstr; strip++) {
	105	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);
	106	hChannel->SetDirectory(0);
	107	sectorArray->AddAtAndExpand(hChannel,strip);
	108	}
	109	}
	110	}
	111	}
	112	}
	113
	114	//_____________________________________________________________________
	115	void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray,
	116	UShort_t det,
	117	Char_t ring,
	118	UShort_t sec,
	119	UShort_t strip) {
	120
	121	TGraphErrors* hChannel = new TGraphErrors();
	122	sectorArray->AddAtAndExpand(hChannel,strip);
	123	}
	124
	125	//_____________________________________________________________________
	126	void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {
	127
	128	UShort_t det = digit->Detector();
	129	Char_t ring = digit->Ring();
	130	UShort_t sec = digit->Sector();
	131	UShort_t strip = digit->Strip();
	132
	133
	134	if(strip%fNumberOfStripsPerChip) return;
	135	Int_t VAchip = strip / fNumberOfStripsPerChip;
	136	TH1S* hChannel = GetChannelHistogram(det,ring,sec,VAchip);
	137	hChannel->Fill(digit->Counts());
	138	UpdatePulseAndADC(det,ring,sec,strip);
	139	}
	140
	141	//_____________________________________________________________________
	142	void AliFMDGainDA::Analyse(UShort_t det,
	143	Char_t ring,
	144	UShort_t sec,
145	UShort_t strip) {
146	TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
147	if(!grChannel->GetN()) {
9c958f2b	148	// AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",det, ring , sec, strip));
f631f26d	149	return;
f631f26d	150	}
9c958f2b	151	TF1 fitFunc("fitFunc","pol1",-10,280);
9c958f2b	152	fitFunc.SetParameters(100,3);
f631f26d	153	grChannel->Fit("fitFunc","Q","Q",0,fHighPulse);
	154	AliFMDParameters* pars = AliFMDParameters::Instance();
	155	UInt_t ddl, board,chip,channel;
	156	pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,channel);
	157
	158	Float_t chi2ndf = 0;
9c958f2b	159	if(fitFunc.GetNDF())
9c958f2b	160	chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
f631f26d	161	ddl = ddl + kBaseDDL;
	162
	163	Int_t relStrip = strip%fNumberOfStripsPerChip;
	164
	165	fOutputFile << ddl << ','
	166	<< board << ','
	167	<< chip << ','
	168	<< channel << ','
	169	<< relStrip << ','
9c958f2b	170	<< fitFunc.GetParameter(1) << ','
9c958f2b	171	<< fitFunc.GetParError(1) << ','
f631f26d	172	<< chi2ndf <<"\n";
	173
	174
	175	if(fSaveHistograms) {
	176
	177	gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,strip));
	178	grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
	179	}
	180
9c958f2b	181
9c958f2b	182
f631f26d	183	}
	184
	185	//_____________________________________________________________________
	186	void AliFMDGainDA::WriteHeaderToFile() {
9c958f2b	187	AliFMDParameters* pars = AliFMDParameters::Instance();
9c958f2b	188	fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
f631f26d	189	fOutputFile.write("# Rcu, Board, Chip, Channel, Strip, Gain, Error, Chi2/NDF \n",59);
	190
	191	}
	192
	193	//_____________________________________________________________________
	194	TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
	195	Char_t ring,
	196	UShort_t sec,
	197	UShort_t strip) {
	198
	199	UShort_t Ring = 1;
	200	if(ring == 'O')
	201	Ring = 0;
	202
	203
	204	TObjArray* detArray = static_cast<TObjArray*>(fGainArray.At(det));
	205	TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
	206	TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
	207	TH1S* hChannel = static_cast<TH1S*>(secArray->At(strip));
	208
	209	return hChannel;
	210	}
	211
	212	//_____________________________________________________________________
	213	TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip) {
	214
	215	UShort_t Ring = 1;
	216	if(ring == 'O')
	217	Ring = 0;
	218
	219
	220	TObjArray* detArray = static_cast<TObjArray*>(fDetectorArray.At(det));
	221	TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
	222	TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
	223	TGraphErrors* hChannel = static_cast<TGraphErrors*>(secArray->At(strip));
	224
	225	return hChannel;
	226	}
	227
	228	//_____________________________________________________________________
	229	void AliFMDGainDA::UpdatePulseAndADC(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip) {
	230
	231	if(strip%fNumberOfStripsPerChip) return;
	232	if(((GetCurrentEvent())%fPulseLength) && GetCurrentEvent()>0) return;
	233
	234	Int_t VAchip = strip/fNumberOfStripsPerChip;
	235	TH1S* hChannel = GetChannelHistogram(det,ring,sec,VAchip);
	236
	237	if(!hChannel->GetEntries()) {
	238	AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",det, ring , sec, strip));
	239	return;
	240	}
	241	Double_t mean = hChannel->GetMean();
	242	Double_t rms = hChannel->GetRMS();
	243	Double_t pulse = (Double_t)fCurrentPulse*fPulseSize;
	244
9c958f2b	245	hChannel->GetXaxis()->SetRangeUser(mean-4rms,mean+4rms);
f631f26d	246
	247	mean = hChannel->GetMean();
	248	rms = hChannel->GetRMS();
	249
	250	Int_t channelNumber = strip + (GetCurrentEvent()-1)/800;
	251
	252	TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
	253
	254	channel->SetPoint(fCurrentPulse,pulse,mean);
	255	channel->SetPointError(fCurrentPulse,0,rms);
	256
	257	if(fSaveHistograms) {
	258	gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,channelNumber));
9c958f2b	259	hChannel->GetXaxis()->SetRange(0,1023);
f631f26d	260	hChannel->Write(Form("hFMD%d%c_%d_%d_pulse_%d",det,ring,sec,channelNumber,fCurrentPulse));
	261	}
	262
	263	hChannel->Reset();
	264
	265	}
	266
	267	//_____________________________________________________________________
	268	void AliFMDGainDA::ResetPulseAndUpdateChannel() {
	269
	270	fCurrentPulse = 0;
	271
	272	}
	273
	274	//_____________________________________________________________________
	275	void AliFMDGainDA::FinishEvent() {
	276
	277	if(GetCurrentEvent()>0 && (GetCurrentEvent()%fPulseLength==0))
	278	fCurrentPulse++;
	279
	280	if(GetCurrentEvent()>0 && (GetCurrentEvent())%fEventsPerChannel==0)
	281	fCurrentPulse = 0;
	282
	283
	284	}
	285	//_____________________________________________________________________
	286	//
	287	//EOF
	288	//