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


ClassImp(AliT0CalibWalk)

//________________________________________________________________
  AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
				      fWalk(0),
				      fAmpLEDRec(0)
{
  //
}

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

}

//________________________________________________________________
AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
  TNamed(calibda),		
  fWalk(0),
  fAmpLEDRec(0)

{
// 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()
{
  //
}
//________________________________________________________________
void AliT0CalibWalk::SetWalk(Int_t ipmt)
{
  //read QTC walk graph from external file

  Int_t mv, ps; 
  Int_t x[70000], y[70000], index[70000];
  Float_t time[10000],amplitude[10000];
  string buffer;
  Bool_t down=false;
  
  const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-Amp.txt");
  ifstream inFile(filename);
  if(!inFile) AliError(Form("Cannot open file %s !",filename));
  
  Int_t i=0;
  while(getline(inFile,buffer)){
    inFile >> ps >> mv;

    x[i]=ps; y[i]=mv;
    i++;
  }
  inFile.close();
  cout<<" number of data "<<i<<endl;
 
  TMath::Sort(i, y, index,down);
  Int_t amp=0, iin=0, isum=0, sum=0;
  Int_t ind=0;
  for (Int_t ii=0; ii<i; ii++)
    {
      ind=index[ii];
      if(y[ind] == amp)
	{
	  sum +=x[ind];
	  iin++;
	}
      else
	{
	  if(iin>0)
	    time[isum] = Float_t (sum/(iin));
	  else
	    time[isum] =Float_t (x[ind]);
	  amplitude[isum] = Float_t (amp);
	  amp=y[ind];
	  iin=0;
	  isum++;
	  sum=0;
	}
    }

  inFile.close();

  TGraph* gr = new TGraph(isum, amplitude, time);
  fWalk.AddAtAndExpand(gr,ipmt);
}

//________________________________________________________________

void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
{
  // read LED walk from external file
  Float_t mv, ps; 
  Float_t x[100], y[100];
  string buffer;
  
  const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-LED.txt");
   ifstream inFile(filename);
  if(!inFile) {AliError(Form("Cannot open file %s !",filename));}
  
  inFile >> mv>>ps;
  Int_t i=0;
  
  while(getline(inFile,buffer)){
    x[i]=mv; y[i]=ps;	
    inFile >> mv >> ps;
    i++;
  }
  inFile.close();
  Float_t y1[100], x1[100];
  for (Int_t ir=0; ir<i; ir++){
    y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
  TGraph* gr = new TGraph(i,y1,x1);
  fAmpLEDRec.AddAtAndExpand(gr,ipmt);
  
}

//________________________________________________________________

void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
{
  //make walk corerction for preprocessor

  TFile *gFile = TFile::Open(laserFile);
  //gSystem->Load("libSpectrum");
 
  Int_t npeaks = 20;
  Int_t sigma=3.;
  Bool_t down = false;

  Int_t index[20];
  Char_t buf1[10], buf2[10],buf3[20];
  
  
  for (Int_t i=0; i<24; i++)
    {
      if(i>11)	{  sprintf(buf1,"T0_A_%i_CFD",i+1-12); }

      if(i<12)	{ sprintf(buf1,"T0_C_%i_CFD",i+1); }

      sprintf(buf2,"CFD_QTC%i",i+1);
      sprintf(buf3,"CFD_LED%i",i+1);
      cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
      TH2F *qtc_cfd = (TH2F*) gFile->Get(buf2);
      TH2F *led_cfd = (TH2F*) gFile->Get(buf3);
      TH1F *cfd = (TH1F*) gFile->Get(buf1);
      //       cfd->Draw();
      TSpectrum *s = new TSpectrum(2*npeaks,1);
      Int_t nfound = s->Search(cfd,sigma," ",0.05);
      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
      if(nfound!=0){
	Float_t *xpeak = s->GetPositionX();
  	TMath::Sort(nfound, xpeak, index,down);
  	Float_t xp = xpeak[index[0]];
	Int_t xbin = cfd->GetXaxis()->FindBin(xp);
	Float_t yp = cfd->GetBinContent(xbin);
	cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
	Float_t hmax = xp+10*sigma;
	Float_t hmin = xp-10*sigma;
	cout<<hmin<< " "<<hmax<<endl;
	//	    cfd->GetXaxis()->SetRange(hmin,hmax);
	//	    TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
	// cfd->Fit("g1","R");
	//	Int_t hminbin=  qtc_cfd->GetXaxis()->GetFirst();
	//	Int_t hmaxbin=  qtc_cfd->GetXaxis()->GetLast();
	Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
	cout<<" nbins "<<nbins<<endl;
	//	cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
	//  qtc_cfd->Draw();
	TProfile *pr_y = qtc_cfd->ProfileX();
	//	Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
	
	pr_y->SetMaximum(hmax);
	pr_y->SetMinimum(hmin);
	Int_t np=nbins/20;
	Double_t *xx = new Double_t[np];
	Double_t *yy = new Double_t[np];
	Int_t ng=0;
	Double_t yg=0;
	for (Int_t ip=1; ip<nbins; ip++)
	  {
	    if(ip%20 != 0 ) {
	      if (pr_y->GetBinContent(ip) !=0)
		yg +=pr_y->GetBinContent(ip);
	      //	cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
	      ng++;}
	    else {
	      xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
	      yy[ip/20] = yg/ng;
	      yg=0;
	      ng=0;
	      //	     	      cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
	    }
	  }
	TGraph *gr = new TGraph(np,xx,yy);
	gr->SetMinimum(hmin);
	gr->SetMaximum(hmax);
	fWalk.AddAtAndExpand(gr,i);	  
	//	fWalk.AddAtAndExpand(gr,i+12);

	Int_t nbinsled= led_cfd->GetXaxis()->GetNbins();
	cout<<" nbins led "<<nbinsled<<endl;
	//	cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
	//  qtc_cfd->Draw();
	TProfile *prled_y = led_cfd->ProfileX();
	//	Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
	
	prled_y->SetMaximum(hmax);
	prled_y->SetMinimum(hmin);
	Int_t npled=nbinsled/20;
	Double_t *xxled = new Double_t[np];
	Double_t *yyled = new Double_t[np];
	Int_t ngled=0;
	Double_t ygled=0;
	for (Int_t ip=1; ip<nbinsled; ip++)
	  {
	    if(ip%20 != 0 ) {
	      if (prled_y->GetBinContent(ip) !=0)
		ygled +=prled_y->GetBinContent(ip);
	      //	cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
	      ngled++;}
	    else {
	      xxled[ip/20] = Float_t (prled_y->GetBinCenter(ip));
	      yyled[ip/20] = ygled/ngled;
	      ygled=0;
	      ngled=0;
	    }
	  }
	TGraph *grled = new TGraph(npled,xxled,yyled);
	grled->SetMinimum(hmin);
	grled->SetMaximum(hmax);
	fAmpLEDRec.AddAtAndExpand(grled,i);	  
 	//	delete [] xx;
	//	delete [] yy;
	
      }
      
    }
  

}
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	**************************************************************************/
	15
	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>
a2ad8166	36
a2ad8166	37
	38	ClassImp(AliT0CalibWalk)
	39
	40	//________________________________________________________________
	41	AliT0CalibWalk::AliT0CalibWalk(): TNamed(),
	42	fWalk(0),
	43	fAmpLEDRec(0)
	44	{
	45	//
	46	}
	47
	48	//________________________________________________________________
	49	AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
	50	fWalk(0),
	51	fAmpLEDRec(0)
	52	{
	53	TString namst = "Calib_";
	54	namst += name;
	55	SetName(namst.Data());
	56	SetTitle(namst.Data());
	57
	58	}
	59
	60	//________________________________________________________________
	61	AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
	62	TNamed(calibda),
	63	fWalk(0),
	64	fAmpLEDRec(0)
	65
	66	{
	67	// copy constructor
	68	SetName(calibda.GetName());
	69	SetTitle(calibda.GetName());
	70
	71
	72	}
	73
	74	//________________________________________________________________
	75	AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
	76	{
	77	// assignment operator
	78	SetName(calibda.GetName());
	79	SetTitle(calibda.GetName());
	80
	81	return *this;
	82	}
	83
	84	//________________________________________________________________
	85	AliT0CalibWalk::~AliT0CalibWalk()
	86	{
	87	//
	88	}
	89	//________________________________________________________________
	90	void AliT0CalibWalk::SetWalk(Int_t ipmt)
	91	{
7576945f	92	//read QTC walk graph from external file
a2ad8166	93
	94	Int_t mv, ps;
	95	Int_t x[70000], y[70000], index[70000];
	96	Float_t time[10000],amplitude[10000];
	97	string buffer;
	98	Bool_t down=false;
	99
	100	const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-Amp.txt");
	101	ifstream inFile(filename);
	102	if(!inFile) AliError(Form("Cannot open file %s !",filename));
	103
	104	Int_t i=0;
	105	while(getline(inFile,buffer)){
	106	inFile >> ps >> mv;
	107
	108	x[i]=ps; y[i]=mv;
	109	i++;
	110	}
	111	inFile.close();
	112	cout<<" number of data "<<i<<endl;
	113
	114	TMath::Sort(i, y, index,down);
	115	Int_t amp=0, iin=0, isum=0, sum=0;
	116	Int_t ind=0;
	117	for (Int_t ii=0; ii<i; ii++)
	118	{
	119	ind=index[ii];
	120	if(y[ind] == amp)
	121	{
	122	sum +=x[ind];
	123	iin++;
	124	}
	125	else
	126	{
	127	if(iin>0)
	128	time[isum] = Float_t (sum/(iin));
	129	else
	130	time[isum] =Float_t (x[ind]);
	131	amplitude[isum] = Float_t (amp);
	132	amp=y[ind];
	133	iin=0;
	134	isum++;
	135	sum=0;
	136	}
	137	}
	138
	139	inFile.close();
	140
	141	TGraph* gr = new TGraph(isum, amplitude, time);
	142	fWalk.AddAtAndExpand(gr,ipmt);
	143	}
	144
	145	//________________________________________________________________
	146
	147	void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
	148	{
7576945f	149	// read LED walk from external file
a2ad8166	150	Float_t mv, ps;
	151	Float_t x[100], y[100];
	152	string buffer;
	153
	154	const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-LED.txt");
	155	ifstream inFile(filename);
	156	if(!inFile) {AliError(Form("Cannot open file %s !",filename));}
	157
	158	inFile >> mv>>ps;
	159	Int_t i=0;
	160
	161	while(getline(inFile,buffer)){
	162	x[i]=mv; y[i]=ps;
	163	inFile >> mv >> ps;
	164	i++;
	165	}
	166	inFile.close();
	167	Float_t y1[100], x1[100];
	168	for (Int_t ir=0; ir<i; ir++){
	169	y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
	170	TGraph* gr = new TGraph(i,y1,x1);
	171	fAmpLEDRec.AddAtAndExpand(gr,ipmt);
	172
	173	}
	174
	175	//________________________________________________________________
	176
	177	void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
	178	{
7576945f	179	//make walk corerction for preprocessor
7576945f	180
a2ad8166	181	TFile *gFile = TFile::Open(laserFile);
a8d3b391	182	//gSystem->Load("libSpectrum");
700e6b36	183
a2ad8166	184	Int_t npeaks = 20;
700e6b36	185	Int_t sigma=3.;
700e6b36	186	Bool_t down = false;
a2ad8166	187
a2ad8166	188	Int_t index[20];
700e6b36	189	Char_t buf1[10], buf2[10],buf3[20];
	190
	191
	192	for (Int_t i=0; i<24; i++)
a2ad8166	193	{
700e6b36	194	if(i>11) { sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
6b7d32f7	195
700e6b36	196	if(i<12) { sprintf(buf1,"T0_C_%i_CFD",i+1); }
f35a2c6c	197
700e6b36	198	sprintf(buf2,"CFD_QTC%i",i+1);
	199	sprintf(buf3,"CFD_LED%i",i+1);
	200	cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
	201	TH2F qtc_cfd = (TH2F) gFile->Get(buf2);
	202	TH2F led_cfd = (TH2F) gFile->Get(buf3);
	203	TH1F cfd = (TH1F) gFile->Get(buf1);
	204	// cfd->Draw();
	205	TSpectrum s = new TSpectrum(2npeaks,1);
	206	Int_t nfound = s->Search(cfd,sigma," ",0.05);
	207	cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
	208	if(nfound!=0){
	209	Float_t *xpeak = s->GetPositionX();
	210	TMath::Sort(nfound, xpeak, index,down);
	211	Float_t xp = xpeak[index[0]];
	212	Int_t xbin = cfd->GetXaxis()->FindBin(xp);
	213	Float_t yp = cfd->GetBinContent(xbin);
	214	cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
	215	Float_t hmax = xp+10*sigma;
	216	Float_t hmin = xp-10*sigma;
	217	cout<<hmin<< " "<<hmax<<endl;
	218	// cfd->GetXaxis()->SetRange(hmin,hmax);
	219	// TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
	220	// cfd->Fit("g1","R");
	221	// Int_t hminbin= qtc_cfd->GetXaxis()->GetFirst();
	222	// Int_t hmaxbin= qtc_cfd->GetXaxis()->GetLast();
	223	Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
	224	cout<<" nbins "<<nbins<<endl;
	225	// cout<<" qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
	226	// qtc_cfd->Draw();
	227	TProfile *pr_y = qtc_cfd->ProfileX();
	228	// Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
	229
	230	pr_y->SetMaximum(hmax);
	231	pr_y->SetMinimum(hmin);
	232	Int_t np=nbins/20;
	233	Double_t *xx = new Double_t[np];
	234	Double_t *yy = new Double_t[np];
	235	Int_t ng=0;
	236	Double_t yg=0;
	237	for (Int_t ip=1; ip<nbins; ip++)
	238	{
	239	if(ip%20 != 0 ) {
	240	if (pr_y->GetBinContent(ip) !=0)
	241	yg +=pr_y->GetBinContent(ip);
	242	// cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
	243	ng++;}
	244	else {
	245	xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
	246	yy[ip/20] = yg/ng;
	247	yg=0;
	248	ng=0;
	249	// cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
	250	}
	251	}
	252	TGraph *gr = new TGraph(np,xx,yy);
	253	gr->SetMinimum(hmin);
	254	gr->SetMaximum(hmax);
	255	fWalk.AddAtAndExpand(gr,i);
	256	// fWalk.AddAtAndExpand(gr,i+12);
	257
	258	Int_t nbinsled= led_cfd->GetXaxis()->GetNbins();
	259	cout<<" nbins led "<<nbinsled<<endl;
	260	// cout<<" qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
	261	// qtc_cfd->Draw();
262	TProfile *prled_y = led_cfd->ProfileX();
263	// Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
264
265	prled_y->SetMaximum(hmax);
266	prled_y->SetMinimum(hmin);
267	Int_t npled=nbinsled/20;
268	Double_t *xxled = new Double_t[np];
269	Double_t *yyled = new Double_t[np];
270	Int_t ngled=0;
271	Double_t ygled=0;
272	for (Int_t ip=1; ip<nbinsled; ip++)
273	{
274	if(ip%20 != 0 ) {
275	if (prled_y->GetBinContent(ip) !=0)
276	ygled +=prled_y->GetBinContent(ip);
277	// cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
278	ngled++;}
279	else {
280	xxled[ip/20] = Float_t (prled_y->GetBinCenter(ip));
281	yyled[ip/20] = ygled/ngled;
282	ygled=0;
283	ngled=0;
284	}
285	}
286	TGraph *grled = new TGraph(npled,xxled,yyled);
287	grled->SetMinimum(hmin);
288	grled->SetMaximum(hmax);
289	fAmpLEDRec.AddAtAndExpand(grled,i);
290	// delete [] xx;
291	// delete [] yy;
292
f35a2c6c	293	}
700e6b36	294
a2ad8166	295	}
700e6b36	296
700e6b36	297
a2ad8166	298
a2ad8166	299	}
a2ad8166	300