[u/mrichter/AliRoot.git] / T0 / AliT0CalibWalk.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
 
/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// class T0 walk correction      Alla Maevskaya alla@inr.ru 20.11.2007        //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include "AliT0CalibWalk.h"
#include "AliLog.h"

#include <TObjArray.h>
#include <TGraph.h>
#include <TFile.h>
#include <TH2F.h> 
#include <TMath.h>
#include <TSystem.h>
#include <Riostream.h>
#include <TSpectrum.h>
#include <TProfile.h>
#include <TF1.h>

using std::cout;
using std::endl;
using std::ifstream;
ClassImp(AliT0CalibWalk)

//________________________________________________________________
  AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
				      fWalk(0),
				      fAmpLEDRec(0),
				      fQTC(0),
                                      fAmpLED(0), 
                                      fCalibByData(kFALSE)
{
  //
}

//________________________________________________________________
AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
				      fWalk(0),
				      fAmpLEDRec(0),				  				      fQTC(0),
				      fAmpLED(0), 
                                      fCalibByData(kFALSE)
    
{
  TString namst = "Calib_";
  namst += name;
  SetName(namst.Data());
  SetTitle(namst.Data());

}

//________________________________________________________________
AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
  TNamed(calibda),		
  fWalk(0),
  fAmpLEDRec(0),
  fQTC(0),
  fAmpLED(0) , 
  fCalibByData(kFALSE)
{
// copy constructor
  SetName(calibda.GetName());
  SetTitle(calibda.GetName());


}

//________________________________________________________________
AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
{
// assignment operator
  SetName(calibda.GetName());
  SetTitle(calibda.GetName());
 
  return *this;
}

//________________________________________________________________
AliT0CalibWalk::~AliT0CalibWalk()
{
  //
}
//________________________________________________________________  

Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
{
  //make walk corerction for preprocessor
  Float_t sigma,cfdmean, qtmean, ledmean;
  Bool_t ok=true;
  Float_t   mips[50];
  cout<<" @@@@ fCalibByData "<<fCalibByData<<endl;

  gFile = TFile::Open(laserFile);
  if(!gFile) {
    AliError("No input laser data found ");
  }
  else
    {
      //      gFile->ls();
      TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
      Int_t nmips=0;
      for (Int_t ibin=0; ibin<2000; ibin++) {
	Float_t bincont = hAmp->GetBinContent(ibin);
	if(bincont>0){ 
	  mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
	  cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
	  nmips++;
	}	
      }    
      /*     
      if (nmips<17) {
	ok=false;
	return ok;
      } 
      */     
      Float_t x1[50], y1[50]; 
      Float_t x2[50], xx2[50],y2[50];
      Float_t xx1[50],yy1[50], xx[50];
      
      Float_t cfd0 = 0;
      
      for (Int_t ii=0; ii<nmips; ii++)
	x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0; 
      
      for (Int_t i=0; i<24; i++)
	{
	  cfd0 = 0;
	  for (Int_t im=0; im<nmips; im++)
	    {	      
	      TString cfd = Form("hCFD%i_%i",i+1,im+1);
	      TString qtc = Form("hQTC%i_%i",i+1,im+1);
	      TString led = Form("hLED%i_%i",i+1,im+1);
	      
	      TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
	      TH1F *hLED = (TH1F*) gFile->Get(led.Data());
	      TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
	      //	      hCFD->SetDirectory(0);
	      //	      hQTC->SetDirectory(0);
	      //	      hLED->SetDirectory(0);
	      if(!hCFD )
 	      	AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
	      if(!hQTC )
 	      	AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
	      if(!hLED)	      
	      	AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
	      if( hCFD && hCFD->GetEntries()<500 ) {
		ok=false;
		printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
		return ok;
	      }
	      if(hCFD && hCFD->GetEntries()>500 ) {
		if( hCFD->GetRMS() >= 1.5) 
		  GetMeanAndSigma(hCFD, cfdmean, sigma);
		else
		  cfdmean = hCFD->GetMean();
		
		Int_t   maxBin = hCFD->GetMaximumBin(); 
		Double_t  meanEstimate = hCFD->GetBinCenter( maxBin); 
		if(TMath::Abs(meanEstimate - cfdmean) > 20 ) cfdmean = meanEstimate; 
		if (im == 0) cfd0 = cfdmean;
		y1[im] =  cfdmean - cfd0;
	      }	
	      if(hQTC && hQTC->GetEntries()>500) {
		GetMeanAndSigma(hQTC, qtmean, sigma);
		
		x1[im] = qtmean;
		if( x1[im] == 0) {
		  ok=false;
		  printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
		  return ok;
		}
	      }
		
	      if( hLED && hLED->GetEntries()>500) {
		GetMeanAndSigma(hLED, ledmean, sigma);
	      }				   
	      else
		{ 
		  //	  ok=false;
		  printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
		  ledmean=0;
		  //  return ok;
		}
	      x2[im] = ledmean;
	      xx2[im] = x2[nmips-im-1];
		
		xx[im]=mips[im];
		if (hQTC) delete  hQTC;
		if (hCFD) delete  hCFD;
		if (hLED) delete  hLED;
		
	    }
	  
	  for (Int_t imi=0; imi<nmips; imi++)
	    {
	      yy1[imi] = Float_t (mips[nmips-imi-1]);
	      xx1[imi] = x2[nmips-imi-1]; 
	      //	      cout<<" LED rev "<<i<<" "<<imi<<" "<<xx1[imi]<<" "<< yy1[imi]<<"nmips-imi-1 = "<< nmips-imi-1<<endl;
	    }
	  
	  if(i==0) cout<<"Making graphs..."<<endl;
	  
      TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
      grwalkqtc->SetTitle(Form("PMT%i",i));
      TGraph *grwalkled = new TGraph (nmips,x2,y1);
      grwalkled->SetTitle(Form("PMT%i",i));
      if(!fCalibByData)
	fWalk.AddAtAndExpand(grwalkqtc,i);
      fAmpLEDRec.AddAtAndExpand(grwalkled,i);
      //	  cout<<" add walk "<<i<<endl;
      
      //fit amplitude graphs to make comparison wth new one	  
      TGraph *grampled = new TGraph (nmips,xx1,yy1);
      TGraph *grqtc = new TGraph (nmips,x1,xx);
      fQTC.AddAtAndExpand(grqtc,i);	 
      fAmpLED.AddAtAndExpand(grampled,i);
      //	  cout<<" add amp "<<i<<endl;
      
      if(i==23)
	cout<<"Graphs created..."<<endl;   
	}
    }
  Float_t xpoint, ypoint, xdata[250], ydata[250];
  Int_t ipmt;
  if(fCalibByData) {
    cout<<" read ingraph "<<endl;
    ifstream ingraph ("calibfit.txt");
    for (Int_t i=0; i<24; i++) 
      {
	for (Int_t ip=0; ip<200; ip++) 
	  {
	    ingraph>>ipmt>>xpoint>>ypoint; 
	    // cout<<i<<" "<<ipmt<<" "<<ip<<" "<< xpoint<<" "<<ypoint<<endl;  
	    xdata[ip]=xpoint;
	    ydata[ip]=ypoint;
	  }
	for (Int_t ip=200; ip<250; ip++) {
	    xdata[ip] =xdata[ip-1]+10;
	    ydata[ip]=ydata[199];
	}
       
	TGraph *grwalkqtc = new TGraph (250,xdata,ydata);
	grwalkqtc->Print();
	grwalkqtc->SetTitle(Form("PMT%i",i));
	fWalk.AddAtAndExpand(grwalkqtc,i);
	for (Int_t ip=0; ip<250; ip++)
	  {
	    xdata[ip]=0;
	    ydata[ip]=0;
	  }
      }
  }

  
return ok;
}

void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma) 
{

  const double window = 2.;  //fit window 
 
  double meanEstimate, sigmaEstimate; 
  int maxBin;
  maxBin        =  hist->GetMaximumBin(); //position of maximum
  meanEstimate  =  hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
  sigmaEstimate = hist->GetRMS();
  TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
  fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
  hist->Fit("fit","RQ","Q");

  mean  = (Float_t) fit->GetParameter(1);
  sigma = (Float_t) fit->GetParameter(2);
  //  printf(" mean %f max %f sigma %f \n",mean, meanEstimate , sigma);  

  delete fit;
}
Commit	Line	Data
a2ad8166	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
271b5bd3	15
a2ad8166	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// class T0 walk correction Alla Maevskaya alla@inr.ru 20.11.2007 //
	21	// //
	22	///////////////////////////////////////////////////////////////////////////////
	23
	24	#include "AliT0CalibWalk.h"
a2ad8166	25	#include "AliLog.h"
a2ad8166	26
	27	#include <TObjArray.h>
	28	#include <TGraph.h>
	29	#include <TFile.h>
a2ad8166	30	#include <TH2F.h>
	31	#include <TMath.h>
	32	#include <TSystem.h>
	33	#include <Riostream.h>
	34	#include <TSpectrum.h>
a2ad8166	35	#include <TProfile.h>
8f1790c4	36	#include <TF1.h>
a2ad8166	37
cdf48647	38	using std::cout;
cdf48647	39	using std::endl;
5e2539d9	40	using std::ifstream;
a2ad8166	41	ClassImp(AliT0CalibWalk)
	42
	43	//________________________________________________________________
	44	AliT0CalibWalk::AliT0CalibWalk(): TNamed(),
	45	fWalk(0),
c883fdf2	46	fAmpLEDRec(0),
c883fdf2	47	fQTC(0),
d0f0bf2a	48	fAmpLED(0),
d0f0bf2a	49	fCalibByData(kFALSE)
a2ad8166	50	{
	51	//
	52	}
	53
	54	//________________________________________________________________
	55	AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
	56	fWalk(0),
c883fdf2	57	fAmpLEDRec(0), fQTC(0),
d0f0bf2a	58	fAmpLED(0),
d0f0bf2a	59	fCalibByData(kFALSE)
c883fdf2	60
a2ad8166	61	{
	62	TString namst = "Calib_";
	63	namst += name;
	64	SetName(namst.Data());
	65	SetTitle(namst.Data());
	66
	67	}
	68
	69	//________________________________________________________________
	70	AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
	71	TNamed(calibda),
	72	fWalk(0),
c883fdf2	73	fAmpLEDRec(0),
c883fdf2	74	fQTC(0),
d0f0bf2a	75	fAmpLED(0) ,
d0f0bf2a	76	fCalibByData(kFALSE)
a2ad8166	77	{
	78	// copy constructor
	79	SetName(calibda.GetName());
	80	SetTitle(calibda.GetName());
	81
	82
	83	}
	84
	85	//________________________________________________________________
	86	AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
	87	{
	88	// assignment operator
	89	SetName(calibda.GetName());
	90	SetTitle(calibda.GetName());
	91
	92	return *this;
	93	}
	94
	95	//________________________________________________________________
	96	AliT0CalibWalk::~AliT0CalibWalk()
	97	{
	98	//
	99	}
d95247c3	100	//________________________________________________________________
a2ad8166	101
d95247c3	102	Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
a2ad8166	103	{
7576945f	104	//make walk corerction for preprocessor
f1799eed	105	Float_t sigma,cfdmean, qtmean, ledmean;
d95247c3	106	Bool_t ok=true;
f1799eed	107	Float_t mips[50];
d0f0bf2a	108	cout<<" @@@@ fCalibByData "<<fCalibByData<<endl;
d0f0bf2a	109
c2337900	110	gFile = TFile::Open(laserFile);
6318dd46	111	if(!gFile) {
	112	AliError("No input laser data found ");
	113	}
	114	else
	115	{
8f1790c4	116	// gFile->ls();
d95247c3	117	TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
d95247c3	118	Int_t nmips=0;
f1799eed	119	for (Int_t ibin=0; ibin<2000; ibin++) {
d95247c3	120	Float_t bincont = hAmp->GetBinContent(ibin);
	121	if(bincont>0){
	122	mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
	123	cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
	124	nmips++;
	125	}
c69fcbc1	126	}
d0f0bf2a	127	/*
f1799eed	128	if (nmips<17) {
07342b9d	129	ok=false;
07342b9d	130	return ok;
f1799eed	131	}
d0f0bf2a	132	*/
c69fcbc1	133	Float_t x1[50], y1[50];
	134	Float_t x2[50], xx2[50],y2[50];
	135	Float_t xx1[50],yy1[50], xx[50];
95f63b36	136
f1799eed	137	Float_t cfd0 = 0;
95f63b36	138
f1799eed	139	for (Int_t ii=0; ii<nmips; ii++)
f1799eed	140	x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0;
95f63b36	141
95f63b36	142	for (Int_t i=0; i<24; i++)
c883fdf2	143	{
f1799eed	144	cfd0 = 0;
f1799eed	145	for (Int_t im=0; im<nmips; im++)
8f1790c4	146	{
d95247c3	147	TString cfd = Form("hCFD%i_%i",i+1,im+1);
d95247c3	148	TString qtc = Form("hQTC%i_%i",i+1,im+1);
3740a735	149	TString led = Form("hLED%i_%i",i+1,im+1);
95f63b36	150
d95247c3	151	TH1F hCFD = (TH1F) gFile->Get(cfd.Data()) ;
	152	TH1F hLED = (TH1F) gFile->Get(led.Data());
	153	TH1F hQTC = (TH1F) gFile->Get(qtc.Data()) ;
f1799eed	154	// hCFD->SetDirectory(0);
	155	// hQTC->SetDirectory(0);
	156	// hLED->SetDirectory(0);
d95247c3	157	if(!hCFD )
	158	AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
	159	if(!hQTC )
	160	AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
	161	if(!hLED)
	162	AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
f1799eed	163	if( hCFD && hCFD->GetEntries()<500 ) {
	164	ok=false;
	165	printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
	166	return ok;
	167	}
	168	if(hCFD && hCFD->GetEntries()>500 ) {
	169	if( hCFD->GetRMS() >= 1.5)
	170	GetMeanAndSigma(hCFD, cfdmean, sigma);
8704b209	171	else
f1799eed	172	cfdmean = hCFD->GetMean();
8704b209	173
f1799eed	174	Int_t maxBin = hCFD->GetMaximumBin();
	175	Double_t meanEstimate = hCFD->GetBinCenter( maxBin);
	176	if(TMath::Abs(meanEstimate - cfdmean) > 20 ) cfdmean = meanEstimate;
	177	if (im == 0) cfd0 = cfdmean;
	178	y1[im] = cfdmean - cfd0;
	179	}
	180	if(hQTC && hQTC->GetEntries()>500) {
	181	GetMeanAndSigma(hQTC, qtmean, sigma);
	182
	183	x1[im] = qtmean;
8704b209	184	if( x1[im] == 0) {
07342b9d	185	ok=false;
8704b209	186	printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
07342b9d	187	return ok;
95f63b36	188	}
c5a24053	189	}
f1799eed	190
	191	if( hLED && hLED->GetEntries()>500) {
	192	GetMeanAndSigma(hLED, ledmean, sigma);
	193	}
	194	else
	195	{
d0f0bf2a	196	// ok=false;
f1799eed	197	printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
d0f0bf2a	198	ledmean=0;
d0f0bf2a	199	// return ok;
f1799eed	200	}
	201	x2[im] = ledmean;
	202	xx2[im] = x2[nmips-im-1];
	203
	204	xx[im]=mips[im];
	205	if (hQTC) delete hQTC;
	206	if (hCFD) delete hCFD;
	207	if (hLED) delete hLED;
	208
e47b556f	209	}
d95247c3	210
8f1790c4	211	for (Int_t imi=0; imi<nmips; imi++)
95f63b36	212	{
8f1790c4	213	yy1[imi] = Float_t (mips[nmips-imi-1]);
8f1790c4	214	xx1[imi] = x2[nmips-imi-1];
d95247c3	215	// cout<<" LED rev "<<i<<" "<<imi<<" "<<xx1[imi]<<" "<< yy1[imi]<<"nmips-imi-1 = "<< nmips-imi-1<<endl;
95f63b36	216	}
95f63b36	217
d95247c3	218	if(i==0) cout<<"Making graphs..."<<endl;
c5a24053	219
f1799eed	220	TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
	221	grwalkqtc->SetTitle(Form("PMT%i",i));
	222	TGraph *grwalkled = new TGraph (nmips,x2,y1);
	223	grwalkled->SetTitle(Form("PMT%i",i));
d0f0bf2a	224	if(!fCalibByData)
d0f0bf2a	225	fWalk.AddAtAndExpand(grwalkqtc,i);
f1799eed	226	fAmpLEDRec.AddAtAndExpand(grwalkled,i);
f1799eed	227	// cout<<" add walk "<<i<<endl;
c5a24053	228
f1799eed	229	//fit amplitude graphs to make comparison wth new one
	230	TGraph *grampled = new TGraph (nmips,xx1,yy1);
	231	TGraph *grqtc = new TGraph (nmips,x1,xx);
	232	fQTC.AddAtAndExpand(grqtc,i);
	233	fAmpLED.AddAtAndExpand(grampled,i);
	234	// cout<<" add amp "<<i<<endl;
	235
	236	if(i==23)
	237	cout<<"Graphs created..."<<endl;
	238	}
d0f0bf2a	239	}
	240	Float_t xpoint, ypoint, xdata[250], ydata[250];
	241	Int_t ipmt;
	242	if(fCalibByData) {
	243	cout<<" read ingraph "<<endl;
	244	ifstream ingraph ("calibfit.txt");
	245	for (Int_t i=0; i<24; i++)
	246	{
	247	for (Int_t ip=0; ip<200; ip++)
	248	{
	249	ingraph>>ipmt>>xpoint>>ypoint;
	250	// cout<<i<<" "<<ipmt<<" "<<ip<<" "<< xpoint<<" "<<ypoint<<endl;
	251	xdata[ip]=xpoint;
	252	ydata[ip]=ypoint;
	253	}
	254	for (Int_t ip=200; ip<250; ip++) {
	255	xdata[ip] =xdata[ip-1]+10;
	256	ydata[ip]=ydata[199];
	257	}
	258
	259	TGraph *grwalkqtc = new TGraph (250,xdata,ydata);
	260	grwalkqtc->Print();
	261	grwalkqtc->SetTitle(Form("PMT%i",i));
	262	fWalk.AddAtAndExpand(grwalkqtc,i);
	263	for (Int_t ip=0; ip<250; ip++)
	264	{
	265	xdata[ip]=0;
	266	ydata[ip]=0;
	267	}
	268	}
	269	}
07342b9d	270
d0f0bf2a	271
d0f0bf2a	272	return ok;
a2ad8166	273	}
f1799eed	274
ff49b4c0	275	void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma)
	276	{
	277
657abde7	278	const double window = 2.; //fit window
ff49b4c0	279
	280	double meanEstimate, sigmaEstimate;
	281	int maxBin;
	282	maxBin = hist->GetMaximumBin(); //position of maximum
	283	meanEstimate = hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
	284	sigmaEstimate = hist->GetRMS();
	285	TF1* fit= new TF1("fit","gaus", meanEstimate - windowsigmaEstimate, meanEstimate + windowsigmaEstimate);
	286	fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
	287	hist->Fit("fit","RQ","Q");
	288
	289	mean = (Float_t) fit->GetParameter(1);
	290	sigma = (Float_t) fit->GetParameter(2);
d0f0bf2a	291	// printf(" mean %f max %f sigma %f \n",mean, meanEstimate , sigma);
ff49b4c0	292
	293	delete fit;
	294	}
c883fdf2	295
	296
	297