[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 "TH1D.h"
#include "TH2D.h"
#include "THnSparse.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(),
    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(),
  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),
  fDebugLevel(calib.fDebugLevel)
{
  //
  // copy constructor
  //
}

AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
  //
  //
  //
  ((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){
  //
  //
  //
  // Thresholds more than 8 tracks with small dip angle
  
  const Int_t kMinLaserTracks = 8;
  const Float_t kThrLaser       = 0.3;
  const Float_t kLaserTgl       = 0.01;

  Int_t ntracks = event->GetNumberOfTracks();
  if (ntracks<kMinLaserTracks) return kFALSE;
  Float_t nlaser=0;
  Float_t nall=0;
  for (Int_t i=0;i<ntracks;++i) {
    AliESDtrack *track=event->GetTrack(i);
    if (!track) continue;
    if (track->GetTPCNcls()<=0) continue; 
    nall++;
    if (TMath::Abs(track->GetTgl())<kLaserTgl) nlaser++;
  }
  if (nlaser>kMinLaserTracks) return kTRUE;
  if (nall>0 && nlaser/nall>kThrLaser) return kTRUE;
  return kFALSE;
}


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){
  //
  // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
  // 
  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){
      TF1 funcGaus("funcGaus","gaus");
      // cut on +- 4 RMS
      projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
      //  
      xvec[counter] = xcenter;
      yvec[counter] = funcGaus.GetParameter(1);
      xerr[counter] = xrms;
      yerr[counter] = funcGaus.GetParError(1); 
      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;
}
Commit	Line	Data
	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	//
	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
	41	//
	42	#include "AliTPCcalibBase.h"
	43	#include "TSystem.h"
	44	#include "TFile.h"
	45	#include "TTreeStream.h"
	46	#include "TTimeStamp.h"
	47	#include "TGraph.h"
	48	#include "TGraphErrors.h"
	49	#include "TF1.h"
	50	#include "TH1D.h"
	51	#include "TH2D.h"
	52	#include "THnSparse.h"
	53
	54
	55	#include "AliLog.h"
	56	#include "AliESDEvent.h"
	57
	58
	59	ClassImp(AliTPCcalibBase)
	60
	61	AliTPCcalibBase::AliTPCcalibBase():
	62	TNamed(),
	63	fDebugStreamer(0),
	64	fStreamLevel(0),
	65	fRun(0), //! current Run number
	66	fEvent(0), //! current Event number
	67	fTime(0), //! current Time
	68	fTrigger(0), //! current trigger type
	69	fMagF(0), //! current magnetic field
	70	fTriggerMaskReject(-1), //trigger mask - reject trigger
	71	fTriggerMaskAccept(-1), //trigger mask - accept trigger
	72	fHasLaser(kFALSE), //flag the laser is overlayed with given event
	73	fRejectLaser(kTRUE), //flag- reject laser
	74	fTriggerClass(),
	75	fDebugLevel(0)
	76	{
	77	//
	78	// Constructor
	79	//
	80	}
	81
	82	AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
	83	TNamed(name,title),
	84	fDebugStreamer(0),
	85	fStreamLevel(0),
	86	fRun(0), //! current Run number
	87	fEvent(0), //! current Event number
	88	fTime(0), //! current Time
	89	fTrigger(0), //! current trigger type
	90	fMagF(0), //! current magnetic field
	91	fTriggerMaskReject(-1), //trigger mask - reject trigger
	92	fTriggerMaskAccept(-1), //trigger mask - accept trigger
	93	fHasLaser(kFALSE), //flag the laser is overlayed with given event
	94	fRejectLaser(kTRUE), //flag- reject laser
	95	fTriggerClass(),
	96	fDebugLevel(0)
	97	{
	98	//
	99	// Constructor
	100	//
	101	}
	102
	103	AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
	104	TNamed(calib),
	105	fDebugStreamer(0),
	106	fStreamLevel(calib.fStreamLevel),
	107	fRun(0), //! current Run number
	108	fEvent(0), //! current Event number
	109	fTime(0), //! current Time
	110	fTrigger(0), //! current trigger type
	111	fMagF(0), //! current magnetic field
	112	fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
	113	fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
	114	fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
	115	fRejectLaser(calib.fRejectLaser), //flag- reject laser
	116	fTriggerClass(calib.fTriggerClass),
	117	fDebugLevel(calib.fDebugLevel)
	118	{
	119	//
	120	// copy constructor
	121	//
	122	}
	123
	124	AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
	125	//
	126	//
	127	//
	128	((TNamed *)this)->operator=(calib);
	129	fDebugStreamer=0;
	130	fStreamLevel=calib.fStreamLevel;
	131	fDebugLevel=calib.fDebugLevel;
	132	return *this;
	133	}
	134
	135
	136	AliTPCcalibBase::~AliTPCcalibBase() {
	137	//
	138	// destructor
	139	//
	140	if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
	141	if (fDebugStreamer) delete fDebugStreamer;
	142	fDebugStreamer=0;
	143	}
	144
	145	void AliTPCcalibBase::Terminate(){
	146	//
	147	//
	148	//
	149	if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
	150	if (fDebugStreamer) delete fDebugStreamer;
	151	fDebugStreamer = 0;
	152	return;
	153	}
	154
	155	TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
	156	//
	157	// Get Debug streamer
	158	// In case debug streamer not yet initialized and StreamLevel>0 create new one
	159	//
	160	if (fStreamLevel==0) return 0;
	161	if (fDebugStreamer) return fDebugStreamer;
	162	TString dsName;
	163	dsName=GetName();
	164	dsName+="Debug.root";
	165	dsName.ReplaceAll(" ","");
	166	fDebugStreamer = new TTreeSRedirector(dsName.Data());
	167	return fDebugStreamer;
	168	}
	169
	170
	171	void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
	172	//
	173	//
	174	//
	175	fRun = event->GetRunNumber();
	176	fEvent = event->GetEventNumberInFile();
	177	fTime = event->GetTimeStamp();
	178	fTrigger = event->GetTriggerMask();
	179	fMagF = event->GetMagneticField();
	180	fTriggerClass = event->GetFiredTriggerClasses().Data();
	181	fHasLaser = HasLaser(event);
	182
	183	}
	184
	185
	186	Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
	187	//
	188	//
	189	//
	190	// Thresholds more than 8 tracks with small dip angle
	191
	192	const Int_t kMinLaserTracks = 8;
	193	const Float_t kThrLaser = 0.3;
	194	const Float_t kLaserTgl = 0.01;
	195
	196	Int_t ntracks = event->GetNumberOfTracks();
	197	if (ntracks<kMinLaserTracks) return kFALSE;
	198	Float_t nlaser=0;
	199	Float_t nall=0;
	200	for (Int_t i=0;i<ntracks;++i) {
	201	AliESDtrack *track=event->GetTrack(i);
	202	if (!track) continue;
	203	if (track->GetTPCNcls()<=0) continue;
	204	nall++;
	205	if (TMath::Abs(track->GetTgl())<kLaserTgl) nlaser++;
	206	}
	207	if (nlaser>kMinLaserTracks) return kTRUE;
	208	if (nall>0 && nlaser/nall>kThrLaser) return kTRUE;
	209	return kFALSE;
	210	}
	211
	212
	213
	214	Bool_t AliTPCcalibBase::AcceptTrigger(){
	215	//
	216	// Apply trigger mask - Don't do calibration for non proper triggers
	217	//
	218	if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
	219	if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
	220	if (fHasLaser && fRejectLaser) return kFALSE;
	221	return kTRUE;
	222	}
	223
	224
	225	void AliTPCcalibBase::RegisterDebugOutput(const char *path){
	226	//
	227	// store - copy debug output to the destination position
	228	// currently ONLY for local copy
	229	if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
	230	if (fStreamLevel==0) return;
	231	TString dsName;
	232	dsName=GetName();
	233	dsName+="Debug.root";
	234	dsName.ReplaceAll(" ","");
	235	TString dsName2=path;
	236	gSystem->MakeDirectory(dsName2.Data());
	237	dsName2+=gSystem->HostName();
	238	gSystem->MakeDirectory(dsName2.Data());
	239	dsName2+="/";
	240	TTimeStamp s;
	241	dsName2+=Int_t(s.GetNanoSec());
	242	dsName2+="/";
	243	gSystem->MakeDirectory(dsName2.Data());
	244	dsName2+=dsName;
	245	AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
	246	printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
	247	TFile::Cp(dsName.Data(),dsName2.Data());
	248	}
	249
	250
	251
	252	TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp){
	253	//
	254	// Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
	255	//
	256	TH2D * hist = h->Projection(axisDim1, axisDim2);
	257	Double_t *xvec = new Double_t[hist->GetNbinsX()];
	258	Double_t *yvec = new Double_t[hist->GetNbinsX()];
	259	Double_t *xerr = new Double_t[hist->GetNbinsX()];
	260	Double_t *yerr = new Double_t[hist->GetNbinsX()];
	261	Int_t counter = 0;
	262	TH1D * projectionHist =0;
	263	//
	264
	265	for(Int_t i=1; i < hist->GetNbinsX(); i++) {
	266	Int_t nsum=0;
	267	Int_t imin = i;
	268	Int_t imax = i;
	269
	270	for (Int_t idelta=0; idelta<nmaxBin; idelta++){
	271	//
	272	imin = TMath::Max(i-idelta,1);
	273	imax = TMath::Min(i+idelta,hist->GetNbinsX());
	274	nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
	275	if (nsum==0) break;
	276	if (nsum>minEntries) break;
	277	}
	278	if (nsum<minEntries) continue;
	279	//
	280	hist->GetXaxis()->SetRange(imin,imax);
	281	projectionHist = hist->ProjectionY("gain",imin,imax);
	282	// Determine Median:
	283	Float_t xMin = projectionHist->GetXaxis()->GetXmin();
	284	Float_t xMax = projectionHist->GetXaxis()->GetXmax();
	285	Float_t xMedian = (xMin+xMax)*0.5;
	286	Float_t integral = 0;
	287	for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
	288	integral+=projectionHist->GetBinContent(jbin);
	289	}
	290	printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
	291	//
	292	//
	293	Float_t currentSum=0;
	294	for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
	295	currentSum += projectionHist->GetBinContent(jbin);
	296	if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
	297	if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
	298	if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
	299	xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
	300	projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
	301	(projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
	302	}
	303	}
	304	//
	305	Float_t rms = projectionHist->GetRMS();
	306	// i += interval;
	307	//
	308	Double_t xcenter = hist->GetMean();
	309	Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
	310	Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
	311	if (rms>0){
	312	TF1 funcGaus("funcGaus","gaus");
	313	// cut on +- 4 RMS
	314	projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
	315	//
	316	xvec[counter] = xcenter;
	317	yvec[counter] = funcGaus.GetParameter(1);
	318	xerr[counter] = xrms;
	319	yerr[counter] = funcGaus.GetParError(1);
	320	counter++;
	321	}else{
	322	xvec[counter] = xcenter;
	323	yvec[counter] = xMedian;
	324	xerr[counter] = xrms;
	325	yerr[counter] = binWidth/TMath::Sqrt(12.);
	326	counter++;
	327	}
	328	delete projectionHist;
	329	}
	330
	331	TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
	332	delete [] xvec;
	333	delete [] yvec;
	334	delete [] xerr;
	335	delete [] yerr;
	336	delete hist;
	337	return graphErrors;
	338	}