[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>

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