[u/mrichter/AliRoot.git] / TPC / AliTPCcalibBase.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.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Base class for the calibration components using 
//  as input TPCseeds and ESDs
//  Event loop outside of the component
//
//
// Base functionality to be implemeneted by component 
/* 
   //In some cases only one of this function to be implemented
   virtual void     Process(AliESDEvent *event)
   virtual void     Process(AliTPCseed *track)
   //
   virtual Long64_t Merge(TCollection *li);
   virtual void     Analyze()
   void             Terminate();
*/
// Functionality provided by base class for Algorith debuging:
//  TTreeSRedirector * cstream =  GetDebugStreamer() - get debug streamer which can be use for numerical debugging
//                      


//  marian.ivanov@cern.ch
// 
#include "AliTPCcalibBase.h"
#include "TSystem.h"
#include "TFile.h"
#include "TTreeStream.h"
#include "TTimeStamp.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TH1.h"
#include "THnSparse.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TAxis.h"
#include "AliRecoParam.h"


#include "AliLog.h"
#include "AliESDEvent.h"


ClassImp(AliTPCcalibBase)

AliTPCcalibBase::AliTPCcalibBase():
    TNamed(),
    fDebugStreamer(0),
    fStreamLevel(0),   
    fRun(0),                  //!  current Run number
    fEvent(0),                //!  current Event number
    fTime(0),                 //!  current Time
    fTrigger(0),              //! current trigger type
    fMagF(0),                 //! current magnetic field
    fTriggerMaskReject(-1),   //trigger mask - reject trigger
    fTriggerMaskAccept(-1),   //trigger mask - accept trigger
    fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
    fRejectLaser(kTRUE),                 //flag- reject laser
    fTriggerClass(),
    fCurrentEvent(0),           //! current event
    fCurrentTrack(0),           //! current esd track
    fCurrentFriendTrack(0),           //! current esd track
    fCurrentSeed(0),            //! current seed
    fDebugLevel(0)
{
  //
  // Constructor
  //
}

AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
  TNamed(name,title),
  fDebugStreamer(0),
  fStreamLevel(0),   
  fRun(0),                  //!  current Run number
  fEvent(0),                //!  current Event number
  fTime(0),                 //!  current Time
  fTrigger(0),              //! current trigger type
  fMagF(0),                 //! current magnetic field
  fTriggerMaskReject(-1),   //trigger mask - reject trigger
  fTriggerMaskAccept(-1),   //trigger mask - accept trigger
  fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
  fRejectLaser(kTRUE),                 //flag- reject laser
  fTriggerClass(),
  fCurrentEvent(0),           //! current event
  fCurrentTrack(0),           //! current esd track
  fCurrentFriendTrack(0),           //! current esd track
  fCurrentSeed(0),            //! current seed
  fDebugLevel(0)
{
  //
  // Constructor
  //
}

AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
  TNamed(calib),
  fDebugStreamer(0),
  fStreamLevel(calib.fStreamLevel),
  fRun(0),                  //!  current Run number
  fEvent(0),                //!  current Event number
  fTime(0),                 //!  current Time
  fTrigger(0),              //! current trigger type
  fMagF(0),                 //! current magnetic field
  fTriggerMaskReject(calib.fTriggerMaskReject),   //trigger mask - reject trigger
  fTriggerMaskAccept(calib.fTriggerMaskAccept),   //trigger mask - accept trigger
  fHasLaser(calib.fHasLaser),                    //flag the laser is overlayed with given event
  fRejectLaser(calib.fRejectLaser),                 //flag- reject laser
  fTriggerClass(calib.fTriggerClass),
  fCurrentEvent(0),           //! current event
  fCurrentTrack(0),           //! current esd track
  fCurrentFriendTrack(0),           //! current esd track
  fCurrentSeed(0),            //! current seed
  fDebugLevel(calib.fDebugLevel)
{
  //
  // copy constructor
  //
}

AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
  //
  // operator=
  //
  if (this == &calib) return (*this);

  ((TNamed *)this)->operator=(calib);
  fDebugStreamer=0;
  fStreamLevel=calib.fStreamLevel;
  fDebugLevel=calib.fDebugLevel;
  return *this;
}


AliTPCcalibBase::~AliTPCcalibBase() {
  //
  // destructor
  //
  if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
  if (fDebugStreamer) delete fDebugStreamer;
  fDebugStreamer=0;
}

void  AliTPCcalibBase::Terminate(){
  //
  //
  //
  if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
  if (fDebugStreamer) delete fDebugStreamer;
  fDebugStreamer = 0;
  return;
}

TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
  //
  // Get Debug streamer
  // In case debug streamer not yet initialized and StreamLevel>0 create new one
  //
  if (fStreamLevel==0) return 0;
  if (fDebugStreamer) return fDebugStreamer;
  TString dsName;
  dsName=GetName();
  dsName+="Debug.root";
  dsName.ReplaceAll(" ",""); 
  fDebugStreamer = new TTreeSRedirector(dsName.Data());
  return fDebugStreamer;
}


void    AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
  //
  //
  //
  fRun     = event->GetRunNumber();
  fEvent   = event->GetEventNumberInFile();
  fTime    = event->GetTimeStamp();
  fTrigger = event->GetTriggerMask();
  fMagF    = event->GetMagneticField();
  fTriggerClass = event->GetFiredTriggerClasses().Data();
  fHasLaser = HasLaser(event); 
  
}


Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
  //
  //
  //
  // Use event specie
  Bool_t isLaser=kFALSE;
  UInt_t specie = event->GetEventSpecie();  // select only cosmic events
  if (specie==AliRecoParam::kCalib) {
    isLaser = kTRUE;
  }
  return isLaser;
}


Bool_t AliTPCcalibBase::AcceptTrigger(){
  //
  // Apply trigger mask - Don't do calibration for non proper triggers
  // 
  if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
  if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
  if (fHasLaser && fRejectLaser) return kFALSE;
  return kTRUE;
}


void AliTPCcalibBase::RegisterDebugOutput(const char *path){
  //
  // store  - copy debug output to the destination position
  // currently ONLY for local copy
  if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
  if (fStreamLevel==0) return;
  TString dsName;
  dsName=GetName();
  dsName+="Debug.root";
  dsName.ReplaceAll(" ",""); 
  TString dsName2=path;
  gSystem->MakeDirectory(dsName2.Data());
  dsName2+=gSystem->HostName();
  gSystem->MakeDirectory(dsName2.Data());
  dsName2+="/";
  TTimeStamp s;
  dsName2+=Int_t(s.GetNanoSec());
  dsName2+="/";
  gSystem->MakeDirectory(dsName2.Data());
  dsName2+=dsName;
  AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
  printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
  TFile::Cp(dsName.Data(),dsName2.Data());
}


TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp, Bool_t useMedian, TTreeSRedirector *cstream, Int_t ival){
  //
  // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
  // 
  TF1 funcGaus("funcGaus","gaus");
  TH2D * hist = h->Projection(axisDim1, axisDim2);
  Double_t *xvec = new Double_t[hist->GetNbinsX()];
  Double_t *yvec = new Double_t[hist->GetNbinsX()];
  Double_t *xerr = new Double_t[hist->GetNbinsX()];
  Double_t *yerr = new Double_t[hist->GetNbinsX()];
  Int_t counter  = 0;
  TH1D * projectionHist =0;
  //

  for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
    Int_t nsum=0;
    Int_t imin   =  i;
    Int_t imax   =  i;

    for (Int_t idelta=0; idelta<nmaxBin; idelta++){
      //
      imin   =  TMath::Max(i-idelta,1);
      imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
      nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
      if (nsum==0) break;
      if (nsum>minEntries) break;
    }
    if (nsum<minEntries) continue;
    //
    hist->GetXaxis()->SetRange(imin,imax);
    projectionHist = hist->ProjectionY("gain",imin,imax);
    // Determine Median:
    Float_t xMin = projectionHist->GetXaxis()->GetXmin();
    Float_t xMax = projectionHist->GetXaxis()->GetXmax();
    Float_t xMedian = (xMin+xMax)*0.5;
    Float_t integral = 0;
    for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
      integral+=projectionHist->GetBinContent(jbin);
    }
    //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
    //
    //
    Float_t currentSum=0;
    for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
      currentSum += projectionHist->GetBinContent(jbin);
      if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
      if (currentSum<fracUp*integral)  xMax = projectionHist->GetBinCenter(jbin+1);      
      if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
	xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
		   projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
	  (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
      }
    }
    //
    Float_t rms  = projectionHist->GetRMS();
    //    i += interval;
    //
    Double_t xcenter =  hist->GetMean(); 
    Double_t xrms    =  hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.); 
    Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
    if (rms>0){
      // cut on +- 4 RMS
      projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
      Double_t chi2 = funcGaus.GetChisquare();
      //  
      xvec[counter] = xcenter;
      yvec[counter] = funcGaus.GetParameter(ival);
      xerr[counter] = xrms;
      yerr[counter] = funcGaus.GetParError(ival); 
      if (useMedian) yvec[counter] = xMedian;
      if (cstream){
	(*cstream)<<"fitDebug"<<
	  "xcenter="<<xcenter<<
	  "xMin="<<xMin<<
	  "xMax="<<xMax<<
	  "xMedian="<<xMedian<<
	  "xFitM"<<yvec[counter]<<
	  "xFitS"<<yerr[counter]<<
	  "chi2="<<chi2<<	  
	  "\n";
      }
      counter++;
    }else{
      xvec[counter] = xcenter;
      yvec[counter] = xMedian;
      xerr[counter] = xrms;
      yerr[counter] = binWidth/TMath::Sqrt(12.); 
      counter++;
    }
    delete projectionHist;
  }
  
  TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
  delete [] xvec;
  delete [] yvec;
  delete [] xerr;
  delete [] yerr;
  delete hist;
  return graphErrors;
}


void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {

  // Method for the correct logarithmic binning of histograms

  TAxis *axis = h->GetAxis(axisDim);
  int bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  Double_t to = axis->GetXmax();
  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = from;
  Double_t factor = pow(to/from, 1./bins);

  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;

}
void AliTPCcalibBase::BinLogX(TH1 *h) {

  // Method for the correct logarithmic binning of histograms

  TAxis *axis = h->GetXaxis();
  int bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  Double_t to = axis->GetXmax();
  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = from;
  Double_t factor = pow(to/from, 1./bins);

  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;

}
Commit	Line	Data
3ab35fc5	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
	17	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// Base class for the calibration components using
	20	// as input TPCseeds and ESDs
	21	// Event loop outside of the component
	22	//
	23	//
ae28e92e	24	// Base functionality to be implemeneted by component
	25	/*
	26	//In some cases only one of this function to be implemented
	27	virtual void Process(AliESDEvent *event)
	28	virtual void Process(AliTPCseed *track)
	29	//
	30	virtual Long64_t Merge(TCollection *li);
	31	virtual void Analyze()
	32	void Terminate();
	33	*/
	34	// Functionality provided by base class for Algorith debuging:
	35	// TTreeSRedirector * cstream = GetDebugStreamer() - get debug streamer which can be use for numerical debugging
	36	//
	37
	38
	39
	40	// marian.ivanov@cern.ch
3ab35fc5	41	//
3ab35fc5	42	#include "AliTPCcalibBase.h"
f7f33dec	43	#include "TSystem.h"
f7f33dec	44	#include "TFile.h"
ae28e92e	45	#include "TTreeStream.h"
47328338	46	#include "TTimeStamp.h"
52a5850f	47	#include "TGraph.h"
	48	#include "TGraphErrors.h"
	49	#include "TF1.h"
3226e48f	50	#include "TH1.h"
3226e48f	51	#include "THnSparse.h"
52a5850f	52	#include "TH1D.h"
52a5850f	53	#include "TH2D.h"
3226e48f	54	#include "TAxis.h"
dcf559ff	55	#include "AliRecoParam.h"
52a5850f	56
	57
	58	#include "AliLog.h"
108953e9	59	#include "AliESDEvent.h"
f7f33dec	60
3ab35fc5	61
	62	ClassImp(AliTPCcalibBase)
	63
ae28e92e	64	AliTPCcalibBase::AliTPCcalibBase():
	65	TNamed(),
	66	fDebugStreamer(0),
108953e9	67	fStreamLevel(0),
	68	fRun(0), //! current Run number
	69	fEvent(0), //! current Event number
	70	fTime(0), //! current Time
	71	fTrigger(0), //! current trigger type
	72	fMagF(0), //! current magnetic field
c51653e8	73	fTriggerMaskReject(-1), //trigger mask - reject trigger
c51653e8	74	fTriggerMaskAccept(-1), //trigger mask - accept trigger
2d22fa36	75	fHasLaser(kFALSE), //flag the laser is overlayed with given event
2d22fa36	76	fRejectLaser(kTRUE), //flag- reject laser
15e48021	77	fTriggerClass(),
b842d904	78	fCurrentEvent(0), //! current event
b842d904	79	fCurrentTrack(0), //! current esd track
4486a91f	80	fCurrentFriendTrack(0), //! current esd track
b842d904	81	fCurrentSeed(0), //! current seed
ae28e92e	82	fDebugLevel(0)
3ab35fc5	83	{
	84	//
	85	// Constructor
	86	//
	87	}
	88
8cf12c96	89	AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
	90	TNamed(name,title),
	91	fDebugStreamer(0),
	92	fStreamLevel(0),
	93	fRun(0), //! current Run number
	94	fEvent(0), //! current Event number
	95	fTime(0), //! current Time
	96	fTrigger(0), //! current trigger type
	97	fMagF(0), //! current magnetic field
	98	fTriggerMaskReject(-1), //trigger mask - reject trigger
	99	fTriggerMaskAccept(-1), //trigger mask - accept trigger
	100	fHasLaser(kFALSE), //flag the laser is overlayed with given event
	101	fRejectLaser(kTRUE), //flag- reject laser
15e48021	102	fTriggerClass(),
b842d904	103	fCurrentEvent(0), //! current event
b842d904	104	fCurrentTrack(0), //! current esd track
4486a91f	105	fCurrentFriendTrack(0), //! current esd track
b842d904	106	fCurrentSeed(0), //! current seed
8cf12c96	107	fDebugLevel(0)
	108	{
	109	//
	110	// Constructor
	111	//
	112	}
	113
ae28e92e	114	AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
	115	TNamed(calib),
	116	fDebugStreamer(0),
	117	fStreamLevel(calib.fStreamLevel),
e5cbd73f	118	fRun(0), //! current Run number
	119	fEvent(0), //! current Event number
	120	fTime(0), //! current Time
	121	fTrigger(0), //! current trigger type
	122	fMagF(0), //! current magnetic field
c51653e8	123	fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
c51653e8	124	fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
2d22fa36	125	fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
2d22fa36	126	fRejectLaser(calib.fRejectLaser), //flag- reject laser
fe8454d9	127	fTriggerClass(calib.fTriggerClass),
b842d904	128	fCurrentEvent(0), //! current event
b842d904	129	fCurrentTrack(0), //! current esd track
4486a91f	130	fCurrentFriendTrack(0), //! current esd track
b842d904	131	fCurrentSeed(0), //! current seed
ae28e92e	132	fDebugLevel(calib.fDebugLevel)
	133	{
	134	//
	135	// copy constructor
	136	//
	137	}
	138
	139	AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
	140	//
4486a91f	141	// operator=
ae28e92e	142	//
04420071	143	if (this == &calib) return (*this);
04420071	144
ae28e92e	145	((TNamed *)this)->operator=(calib);
	146	fDebugStreamer=0;
	147	fStreamLevel=calib.fStreamLevel;
	148	fDebugLevel=calib.fDebugLevel;
525062ee	149	return *this;
ae28e92e	150	}
	151
	152
3ab35fc5	153	AliTPCcalibBase::~AliTPCcalibBase() {
	154	//
	155	// destructor
	156	//
f7f33dec	157	if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
ae28e92e	158	if (fDebugStreamer) delete fDebugStreamer;
	159	fDebugStreamer=0;
	160	}
	161
	162	void AliTPCcalibBase::Terminate(){
	163	//
	164	//
	165	//
f7f33dec	166	if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
ae28e92e	167	if (fDebugStreamer) delete fDebugStreamer;
	168	fDebugStreamer = 0;
	169	return;
	170	}
	171
	172	TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
	173	//
	174	// Get Debug streamer
	175	// In case debug streamer not yet initialized and StreamLevel>0 create new one
	176	//
	177	if (fStreamLevel==0) return 0;
	178	if (fDebugStreamer) return fDebugStreamer;
	179	TString dsName;
	180	dsName=GetName();
	181	dsName+="Debug.root";
	182	dsName.ReplaceAll(" ","");
	183	fDebugStreamer = new TTreeSRedirector(dsName.Data());
	184	return fDebugStreamer;
3ab35fc5	185	}
f7f33dec	186
f7f33dec	187
108953e9	188	void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
	189	//
	190	//
	191	//
	192	fRun = event->GetRunNumber();
	193	fEvent = event->GetEventNumberInFile();
	194	fTime = event->GetTimeStamp();
	195	fTrigger = event->GetTriggerMask();
	196	fMagF = event->GetMagneticField();
15e48021	197	fTriggerClass = event->GetFiredTriggerClasses().Data();
2d22fa36	198	fHasLaser = HasLaser(event);
15e48021	199
108953e9	200	}
108953e9	201
2d22fa36	202
	203	Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
	204	//
	205	//
	206	//
dcf559ff	207	// Use event specie
	208	Bool_t isLaser=kFALSE;
	209	UInt_t specie = event->GetEventSpecie(); // select only cosmic events
	210	if (specie==AliRecoParam::kCalib) {
	211	isLaser = kTRUE;
2d22fa36	212	}
dcf559ff	213	return isLaser;
2d22fa36	214	}
	215
	216
	217
c51653e8	218	Bool_t AliTPCcalibBase::AcceptTrigger(){
	219	//
	220	// Apply trigger mask - Don't do calibration for non proper triggers
	221	//
2d22fa36	222	if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
	223	if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
	224	if (fHasLaser && fRejectLaser) return kFALSE;
c51653e8	225	return kTRUE;
	226	}
	227
	228
f7f33dec	229	void AliTPCcalibBase::RegisterDebugOutput(const char *path){
	230	//
	231	// store - copy debug output to the destination position
	232	// currently ONLY for local copy
	233	if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
	234	if (fStreamLevel==0) return;
	235	TString dsName;
	236	dsName=GetName();
	237	dsName+="Debug.root";
	238	dsName.ReplaceAll(" ","");
	239	TString dsName2=path;
	240	gSystem->MakeDirectory(dsName2.Data());
	241	dsName2+=gSystem->HostName();
	242	gSystem->MakeDirectory(dsName2.Data());
	243	dsName2+="/";
47328338	244	TTimeStamp s;
	245	dsName2+=Int_t(s.GetNanoSec());
	246	dsName2+="/";
	247	gSystem->MakeDirectory(dsName2.Data());
f7f33dec	248	dsName2+=dsName;
	249	AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
	250	printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
	251	TFile::Cp(dsName.Data(),dsName2.Data());
	252	}
52a5850f	253
	254
	255
70b0f0bf	256	TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp, Bool_t useMedian, TTreeSRedirector cstream, Int_t ival){
52a5850f	257	//
	258	// Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
	259	//
3d549531	260	TF1 funcGaus("funcGaus","gaus");
52a5850f	261	TH2D * hist = h->Projection(axisDim1, axisDim2);
	262	Double_t *xvec = new Double_t[hist->GetNbinsX()];
	263	Double_t *yvec = new Double_t[hist->GetNbinsX()];
	264	Double_t *xerr = new Double_t[hist->GetNbinsX()];
	265	Double_t *yerr = new Double_t[hist->GetNbinsX()];
	266	Int_t counter = 0;
	267	TH1D * projectionHist =0;
	268	//
	269
6b8ef6d1	270	for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
52a5850f	271	Int_t nsum=0;
	272	Int_t imin = i;
	273	Int_t imax = i;
	274
	275	for (Int_t idelta=0; idelta<nmaxBin; idelta++){
	276	//
	277	imin = TMath::Max(i-idelta,1);
	278	imax = TMath::Min(i+idelta,hist->GetNbinsX());
	279	nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
	280	if (nsum==0) break;
	281	if (nsum>minEntries) break;
	282	}
	283	if (nsum<minEntries) continue;
	284	//
	285	hist->GetXaxis()->SetRange(imin,imax);
	286	projectionHist = hist->ProjectionY("gain",imin,imax);
	287	// Determine Median:
	288	Float_t xMin = projectionHist->GetXaxis()->GetXmin();
	289	Float_t xMax = projectionHist->GetXaxis()->GetXmax();
	290	Float_t xMedian = (xMin+xMax)*0.5;
d0b31f57	291	Float_t integral = 0;
	292	for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
	293	integral+=projectionHist->GetBinContent(jbin);
	294	}
7f766801	295	//printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
52a5850f	296	//
	297	//
	298	Float_t currentSum=0;
	299	for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
	300	currentSum += projectionHist->GetBinContent(jbin);
d0b31f57	301	if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
	302	if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
	303	if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
52a5850f	304	xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
	305	projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
	306	(projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
	307	}
	308	}
	309	//
	310	Float_t rms = projectionHist->GetRMS();
	311	// i += interval;
	312	//
	313	Double_t xcenter = hist->GetMean();
	314	Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
	315	Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
	316	if (rms>0){
52a5850f	317	// cut on +- 4 RMS
52a5850f	318	projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
3d549531	319	Double_t chi2 = funcGaus.GetChisquare();
52a5850f	320	//
52a5850f	321	xvec[counter] = xcenter;
70b0f0bf	322	yvec[counter] = funcGaus.GetParameter(ival);
52a5850f	323	xerr[counter] = xrms;
70b0f0bf	324	yerr[counter] = funcGaus.GetParError(ival);
3d549531	325	if (useMedian) yvec[counter] = xMedian;
	326	if (cstream){
	327	(*cstream)<<"fitDebug"<<
	328	"xcenter="<<xcenter<<
	329	"xMin="<<xMin<<
	330	"xMax="<<xMax<<
	331	"xMedian="<<xMedian<<
	332	"xFitM"<<yvec[counter]<<
	333	"xFitS"<<yerr[counter]<<
	334	"chi2="<<chi2<<
	335	"\n";
	336	}
52a5850f	337	counter++;
	338	}else{
	339	xvec[counter] = xcenter;
	340	yvec[counter] = xMedian;
	341	xerr[counter] = xrms;
	342	yerr[counter] = binWidth/TMath::Sqrt(12.);
	343	counter++;
	344	}
	345	delete projectionHist;
	346	}
	347
	348	TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
	349	delete [] xvec;
	350	delete [] yvec;
	351	delete [] xerr;
	352	delete [] yerr;
	353	delete hist;
	354	return graphErrors;
	355	}
3226e48f	356
	357
	358	void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
	359
	360	// Method for the correct logarithmic binning of histograms
	361
	362	TAxis *axis = h->GetAxis(axisDim);
	363	int bins = axis->GetNbins();
	364
	365	Double_t from = axis->GetXmin();
	366	Double_t to = axis->GetXmax();
	367	Double_t *new_bins = new Double_t[bins + 1];
	368
	369	new_bins[0] = from;
	370	Double_t factor = pow(to/from, 1./bins);
	371
	372	for (int i = 1; i <= bins; i++) {
	373	new_bins[i] = factor * new_bins[i-1];
	374	}
	375	axis->Set(bins, new_bins);
4ce766eb	376	delete [] new_bins;
3226e48f	377
	378	}
	379	void AliTPCcalibBase::BinLogX(TH1 *h) {
	380
	381	// Method for the correct logarithmic binning of histograms
	382
	383	TAxis *axis = h->GetXaxis();
	384	int bins = axis->GetNbins();
	385
	386	Double_t from = axis->GetXmin();
	387	Double_t to = axis->GetXmax();
	388	Double_t *new_bins = new Double_t[bins + 1];
	389
	390	new_bins[0] = from;
	391	Double_t factor = pow(to/from, 1./bins);
	392
	393	for (int i = 1; i <= bins; i++) {
	394	new_bins[i] = factor * new_bins[i-1];
	395	}
	396	axis->Set(bins, new_bins);
4ce766eb	397	delete [] new_bins;
3226e48f	398
3226e48f	399	}