[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDpidRefMakerLQ.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: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */

///////////////////////////////////////////////////////////////////////////////
//
//
//  TRD calibration class for building reference data for PID
//  - 2D reference histograms (responsible A.Bercuci) 
//  - 3D reference histograms (not yet implemented) (responsible A.Bercuci)
//
//   Origin
//   Alex Bercuci  (A.Bercuci@gsi.de)
//
///////////////////////////////////////////////////////////////////////////////

#include <TSystem.h>
#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TMath.h>
#include <TEventList.h>
#include <TH2D.h>
#include <TH2I.h>
#include <TCanvas.h>

#include "AliLog.h"
#include "../STAT/TKDPDF.h"
#include "AliTRDpidRefMakerLQ.h"
#include "Cal/AliTRDCalPID.h"
#include "Cal/AliTRDCalPIDLQ.h"
#include "AliTRDseedV1.h"
#include "AliTRDcalibDB.h"
#include "AliTRDgeometry.h"

ClassImp(AliTRDpidRefMakerLQ)

//__________________________________________________________________
AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() :
  AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
{
  //
  // AliTRDpidRefMakerLQ default constructor
  //
}

//__________________________________________________________________
AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
{
  //
  // AliTRDCalPIDQRef destructor
  //
}

// //________________________________________________________________________
void AliTRDpidRefMakerLQ::CreateOutputObjects()
{
  // Create histograms
  // Called once

  //OpenFile(0, "RECREATE");
  fContainer = Histos();
}


//__________________________________________________________________
TObjArray* AliTRDpidRefMakerLQ::Histos()
{
  // Create histograms

  if(fContainer) return fContainer;

  fContainer = new TObjArray(AliTRDCalPID::kNMom);
  //fContainer->SetOwner(kTRUE);
  //fContainer->SetName(Form("Moni%s", GetName()));

  // save dE/dx references
  TH2 *h2 = 0x0;
  for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ 
    TObjArray *arr = new TObjArray(AliPID::kSPECIES);
    arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
    for(Int_t is=AliPID::kSPECIES; is--;) {
      h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
      h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
      h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
      h2->GetZaxis()->SetTitle("#");
      arr->AddAt(h2, is);
    }
    fContainer->AddAt(arr, ip);
  }
  fNRefFigures=AliTRDCalPID::kNMom;
  return fContainer;
}


//__________________________________________________________________
TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
{
// Steer loading of OCDB LQ PID

  TDirectoryFile *d = 0x0;
  if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
  if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
  AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
  cal->LoadPDF(d);
  return cal;
}

//__________________________________________________________________
Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
{
// Steering reference picture

  if(ifig<0 || ifig>=fNRefFigures){ 
    AliError("Ref fig requested outside definition.");
    return kFALSE;
  }
  if(!gPad){
    AliWarning("Please provide a canvas to draw results.");
    return kFALSE;
  }

  TObjArray *arr(NULL);TList *l(NULL);TH2 *h2(NULL);
  switch(ifig){
  case 0: // PDG plot
    if(!(h2=(TH2*)fContainer->At(0))){
      AliError("Abundance Plot missing.");
      return kFALSE;
    }
    h2->Draw("cont4z");
    break;
  default:
    if(!(arr = (TObjArray*)fContainer->At(ifig))){
      AliError(Form("2D container @ pBin[%d] missing.", ifig-1));
      return kFALSE;
    }
    gPad->Divide(3, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives(); 
    for(Int_t is=0; is<AliPID::kSPECIES; is++){
      ((TVirtualPad*)l->At(is))->cd();
      if(!(h2=(TH2*)arr->At(is))){
        AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig-1));
        return kFALSE;
      }
      h2->Draw("cont4z");
    }
  }
  return kTRUE;
}


//________________________________________________________________________
Bool_t AliTRDpidRefMakerLQ::Load(const Char_t */*fname*/)
{
  const Char_t *name("PIDrefMaker");
  if(gSystem->AccessPathName(Form("TRD.Calib%s.root", name), kReadPermission)){
    AliError(Form("File TRD.Calib%s.root not readable", name));
    return kFALSE;
  }
  if(!TFile::Open(Form("TRD.Calib%s.root", name))){
    AliError(Form("File TRD.Calib%s.root corrupted", name));
    return kFALSE;
  }
  if (!(fData = dynamic_cast<TTree*>(gFile->Get(name)))) {
    AliError(Form("Tree %s not available", name));
    return kFALSE;
  }
  LinkPIDdata();

  TObjArray *o(NULL);
  if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", name)))){
    AliWarning(Form("Monitor container Moni%s not available.", name));
    return kFALSE;
  }
  fContainer = (TObjArray*)o->Clone("monitor");
  return kTRUE;
}

//________________________________________________________________________
void AliTRDpidRefMakerLQ::Exec(Option_t */*opt*/)
{
// Mock up function to load PID data into local data storage
  AliInfo(Form("fInfo[%d]\n", fInfo->GetEntriesFast()));
  PostData(0, fContainer);
}


//________________________________________________________________________
Bool_t AliTRDpidRefMakerLQ::PostProcess()
{
// Analyse merged dedx = f(p) distributions.
//   - select momentum - species bins
//   - rotate to principal components
//   - locally interpolate with TKDPDF
//   - save interpolation to monitoring histograms
//   - write pdf to file for loading to OCDB
// 

  TDatime d;
  TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
  pdfs->Write();
  AliDebug(2, Form("Data[%d]", fData->GetEntries()));
  pdfs->cd();

  //TCanvas *cc = new TCanvas("cc", "", 500, 500);
  // allocate working storage
  Float_t *data[] = {
    new Float_t[kMaxStat], 
    new Float_t[kMaxStat]};
  for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){ 
    for(Int_t is=AliPID::kSPECIES; is--;) {
      Int_t n(0); // index of data
      for(Int_t itrk=0; (itrk < fData->GetEntries()) && (n<kMaxStat); itrk++){
        if(!(fData->GetEntry(itrk))) continue;
        if(fPIDbin!=is) continue;
        for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
          if((fPIDdataArray->fData[ily].fPLbin & 0xf)!= ip) continue;
          
          Float_t dedx[] = {0., 0.};
          for(Int_t islice=AliTRDCalPID::kNSlicesNN; islice--;){
            Int_t jslice = islice>kNN2LQtransition;
            dedx[jslice]+=fPIDdataArray->fData[ily].fdEdx[islice];
          }
          
          // check data integrity
          if(dedx[0]<1.e-30) continue;
          if(dedx[1]<1.e-30) continue;

          // store data
          data[0][n] = TMath::Log(dedx[0]);
          data[1][n] = TMath::Log(dedx[1]);
          n++; if(n==kMaxStat) break;
        }
      }

      // estimate bucket statistics
      Int_t nb(kMinBuckets), // number of buckets
            ns(Int_t(Float_t(n)/nb));    //statistics/bucket
            
// if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
//       else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error

      AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, n, ns, nb));
      if(ns<Int_t(kMinStat)){
        AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
        continue;
      }

      // build PDF
      TKDPDF pdf(n, 2, ns, data);
      pdf.SetCOG(kFALSE);
      pdf.SetWeights();
      pdf.SetStore();
      pdf.SetAlpha(5.);
      pdf.GetStatus();
      Float_t *c, v, ve; Double_t r, e, rxy[2];
      for(Int_t in=pdf.GetNTNodes(); in--;){
        pdf.GetCOGPoint(in, c, v, ve);
        rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
        pdf.Eval(rxy, r, e, kTRUE);
      }
//       // visual on-line monitoring
//       pdf.DrawProjection();cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
//       cc->SaveAs(Form("%s_%s%02d.gif", GetName(), AliPID::ParticleShortName(is), ip));

      // save a discretization of the PDF for result monitoring
      TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(is);
      TAxis *ax = h2s->GetXaxis(), *ay = h2s->GetYaxis();
      h2s->Clear();
      for(int ix=1; ix<=ax->GetNbins(); ix++){
        rxy[0] = ax->GetBinCenter(ix);
        for(int iy=1; iy<=ay->GetNbins(); iy++){
          rxy[1] = ay->GetBinCenter(iy);
      
          Double_t rr,ee;
          pdf.Eval(rxy, rr, ee, kFALSE);
          if(rr<0. || ee/rr>.15) continue; // 15% relative error
          //printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
          h2s->SetBinContent(ix, iy, rr);
        }
      }

      // write results to output array
      //pdf.GetStatus();
      pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
    }
  }
  delete [] data[0]; delete [] data[1];
  pdfs->Write();
  //fCalib->Close(); delete fCalib;

  return kTRUE; // testing protection
}
Commit	Line	Data
1ee39b3a	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: AliTRDpidRefMakerLQ.cxx 34163 2009-08-07 11:28:51Z cblume $ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	//
	20	//
	21	// TRD calibration class for building reference data for PID
	22	// - 2D reference histograms (responsible A.Bercuci)
	23	// - 3D reference histograms (not yet implemented) (responsible A.Bercuci)
1ee39b3a	24	//
	25	// Origin
	26	// Alex Bercuci (A.Bercuci@gsi.de)
	27	//
	28	///////////////////////////////////////////////////////////////////////////////
	29
b91fdd71	30	#include <TSystem.h>
1ee39b3a	31	#include <TROOT.h>
b91fdd71	32	#include <TFile.h>
1ee39b3a	33	#include <TTree.h>
	34	#include <TMath.h>
	35	#include <TEventList.h>
	36	#include <TH2D.h>
	37	#include <TH2I.h>
1ee39b3a	38	#include <TCanvas.h>
1ee39b3a	39
1ee39b3a	40	#include "AliLog.h"
da27d983	41	#include "../STAT/TKDPDF.h"
1ee39b3a	42	#include "AliTRDpidRefMakerLQ.h"
da27d983	43	#include "Cal/AliTRDCalPID.h"
da27d983	44	#include "Cal/AliTRDCalPIDLQ.h"
1ee39b3a	45	#include "AliTRDseedV1.h"
	46	#include "AliTRDcalibDB.h"
	47	#include "AliTRDgeometry.h"
	48
	49	ClassImp(AliTRDpidRefMakerLQ)
	50
	51	//__________________________________________________________________
b91fdd71	52	AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() :
b91fdd71	53	AliTRDpidRefMaker("PIDrefMakerLQ", "PID(LQ) Reference Maker")
1ee39b3a	54	{
	55	//
	56	// AliTRDpidRefMakerLQ default constructor
	57	//
	58	}
	59
	60	//__________________________________________________________________
	61	AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
	62	{
	63	//
	64	// AliTRDCalPIDQRef destructor
	65	//
	66	}
	67
b91fdd71	68	// //________________________________________________________________________
1ee39b3a	69	void AliTRDpidRefMakerLQ::CreateOutputObjects()
	70	{
	71	// Create histograms
	72	// Called once
	73
b91fdd71	74	//OpenFile(0, "RECREATE");
	75	fContainer = Histos();
	76	}
	77
	78
	79	//__________________________________________________________________
	80	TObjArray* AliTRDpidRefMakerLQ::Histos()
	81	{
	82	// Create histograms
	83
	84	if(fContainer) return fContainer;
	85
	86	fContainer = new TObjArray(AliTRDCalPID::kNMom);
	87	//fContainer->SetOwner(kTRUE);
	88	//fContainer->SetName(Form("Moni%s", GetName()));
1ee39b3a	89
	90	// save dE/dx references
	91	TH2 *h2 = 0x0;
	92	for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
	93	TObjArray *arr = new TObjArray(AliPID::kSPECIES);
b91fdd71	94	arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
1ee39b3a	95	for(Int_t is=AliPID::kSPECIES; is--;) {
b91fdd71	96	h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 5., 10., 50, 5., 10.);
1ee39b3a	97	h2->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
	98	h2->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
	99	h2->GetZaxis()->SetTitle("#");
	100	arr->AddAt(h2, is);
	101	}
b91fdd71	102	fContainer->AddAt(arr, ip);
1ee39b3a	103	}
b91fdd71	104	fNRefFigures=AliTRDCalPID::kNMom;
b91fdd71	105	return fContainer;
1ee39b3a	106	}
1ee39b3a	107
1ee39b3a	108
	109	//__________________________________________________________________
	110	TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t *opt)
	111	{
	112	// Steer loading of OCDB LQ PID
	113
	114	TDirectoryFile *d = 0x0;
	115	if(!TFile::Open(Form("TRD.Calib%s.root", GetName()))) return 0x0;
	116	if(!(d=(TDirectoryFile*)gFile->Get(Form("PDF_%s", opt)))) return 0x0;
	117	AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
	118	cal->LoadPDF(d);
	119	return cal;
	120	}
	121
	122	//__________________________________________________________________
	123	Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
	124	{
	125	// Steering reference picture
	126
b91fdd71	127	if(ifig<0 \|\| ifig>=fNRefFigures){
1ee39b3a	128	AliError("Ref fig requested outside definition.");
	129	return kFALSE;
	130	}
	131	if(!gPad){
	132	AliWarning("Please provide a canvas to draw results.");
	133	return kFALSE;
	134	}
	135
b91fdd71	136	TObjArray arr(NULL);TList l(NULL);TH2 *h2(NULL);
	137	switch(ifig){
	138	case 0: // PDG plot
	139	if(!(h2=(TH2*)fContainer->At(0))){
	140	AliError("Abundance Plot missing.");
	141	return kFALSE;
	142	}
	143	h2->Draw("cont4z");
	144	break;
	145	default:
	146	if(!(arr = (TObjArray*)fContainer->At(ifig))){
	147	AliError(Form("2D container @ pBin[%d] missing.", ifig-1));
	148	return kFALSE;
	149	}
	150	gPad->Divide(3, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
	151	for(Int_t is=0; is<AliPID::kSPECIES; is++){
	152	((TVirtualPad*)l->At(is))->cd();
	153	if(!(h2=(TH2*)arr->At(is))){
	154	AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig-1));
	155	return kFALSE;
	156	}
	157	h2->Draw("cont4z");
	158	}
	159	}
	160	return kTRUE;
	161	}
	162
	163
	164	//________________________________________________________________________
	165	Bool_t AliTRDpidRefMakerLQ::Load(const Char_t /fname*/)
	166	{
	167	const Char_t *name("PIDrefMaker");
	168	if(gSystem->AccessPathName(Form("TRD.Calib%s.root", name), kReadPermission)){
	169	AliError(Form("File TRD.Calib%s.root not readable", name));
	170	return kFALSE;
	171	}
	172	if(!TFile::Open(Form("TRD.Calib%s.root", name))){
	173	AliError(Form("File TRD.Calib%s.root corrupted", name));
	174	return kFALSE;
	175	}
	176	if (!(fData = dynamic_cast<TTree*>(gFile->Get(name)))) {
	177	AliError(Form("Tree %s not available", name));
	178	return kFALSE;
1ee39b3a	179	}
b91fdd71	180	LinkPIDdata();
1ee39b3a	181
b91fdd71	182	TObjArray *o(NULL);
	183	if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", name)))){
	184	AliWarning(Form("Monitor container Moni%s not available.", name));
	185	return kFALSE;
	186	}
	187	fContainer = (TObjArray*)o->Clone("monitor");
1ee39b3a	188	return kTRUE;
	189	}
	190
	191	//________________________________________________________________________
b91fdd71	192	void AliTRDpidRefMakerLQ::Exec(Option_t /opt*/)
	193	{
	194	// Mock up function to load PID data into local data storage
	195	AliInfo(Form("fInfo[%d]\n", fInfo->GetEntriesFast()));
	196	PostData(0, fContainer);
	197	}
	198
	199
	200	//________________________________________________________________________
1ee39b3a	201	Bool_t AliTRDpidRefMakerLQ::PostProcess()
	202	{
	203	// Analyse merged dedx = f(p) distributions.
	204	// - select momentum - species bins
	205	// - rotate to principal components
	206	// - locally interpolate with TKDPDF
	207	// - save interpolation to monitoring histograms
	208	// - write pdf to file for loading to OCDB
	209	//
	210
1ee39b3a	211	TDatime d;
	212	TDirectoryFile *pdfs = new TDirectoryFile(Form("PDF_%d", d.GetDate()), "PDFs for LQ TRD-PID", "", gFile);
	213	pdfs->Write();
	214	AliDebug(2, Form("Data[%d]", fData->GetEntries()));
	215	pdfs->cd();
	216
	217	//TCanvas *cc = new TCanvas("cc", "", 500, 500);
b91fdd71	218	// allocate working storage
	219	Float_t *data[] = {
	220	new Float_t[kMaxStat],
	221	new Float_t[kMaxStat]};
1ee39b3a	222	for(Int_t ip=AliTRDCalPID::kNMom; ip--; ){
	223	for(Int_t is=AliPID::kSPECIES; is--;) {
	224	Int_t n(0); // index of data
	225	for(Int_t itrk=0; (itrk < fData->GetEntries()) && (n<kMaxStat); itrk++){
	226	if(!(fData->GetEntry(itrk))) continue;
	227	if(fPIDbin!=is) continue;
	228	for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
	229	if((fPIDdataArray->fData[ily].fPLbin & 0xf)!= ip) continue;
	230
	231	Float_t dedx[] = {0., 0.};
	232	for(Int_t islice=AliTRDCalPID::kNSlicesNN; islice--;){
	233	Int_t jslice = islice>kNN2LQtransition;
	234	dedx[jslice]+=fPIDdataArray->fData[ily].fdEdx[islice];
	235	}
	236
	237	// check data integrity
	238	if(dedx[0]<1.e-30) continue;
	239	if(dedx[1]<1.e-30) continue;
	240
	241	// store data
	242	data[0][n] = TMath::Log(dedx[0]);
	243	data[1][n] = TMath::Log(dedx[1]);
	244	n++; if(n==kMaxStat) break;
	245	}
	246	}
	247
	248	// estimate bucket statistics
	249	Int_t nb(kMinBuckets), // number of buckets
	250	ns(Int_t(Float_t(n)/nb)); //statistics/bucket
	251
	252	// if(Float_t(n)/nb < 220.) ns = 200; // 7% stat error
	253	// else if(Float_t(n)/nb < 420.) ns = 400; // 5% stat error
	254
	255	AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, n, ns, nb));
	256	if(ns<Int_t(kMinStat)){
	257	AliWarning(Form("Not enough entries [%d] for %s[%d].", n, AliPID::ParticleShortName(is), ip));
	258	continue;
	259	}
	260
	261	// build PDF
	262	TKDPDF pdf(n, 2, ns, data);
	263	pdf.SetCOG(kFALSE);
	264	pdf.SetWeights();
	265	pdf.SetStore();
	266	pdf.SetAlpha(5.);
	267	pdf.GetStatus();
	268	Float_t *c, v, ve; Double_t r, e, rxy[2];
	269	for(Int_t in=pdf.GetNTNodes(); in--;){
	270	pdf.GetCOGPoint(in, c, v, ve);
	271	rxy[0] = (Double_t)c[0];rxy[1] = (Double_t)c[1];
	272	pdf.Eval(rxy, r, e, kTRUE);
	273	}
	274	// // visual on-line monitoring
	275	// pdf.DrawProjection();cc->Modified(); cc->Update(); cc->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
	276	// cc->SaveAs(Form("%s_%s%02d.gif", GetName(), AliPID::ParticleShortName(is), ip));
	277
	278	// save a discretization of the PDF for result monitoring
	279	TH2 h2s = (TH2D)((TObjArray*)fContainer->At(ip))->At(is);
	280	TAxis ax = h2s->GetXaxis(), ay = h2s->GetYaxis();
	281	h2s->Clear();
	282	for(int ix=1; ix<=ax->GetNbins(); ix++){
	283	rxy[0] = ax->GetBinCenter(ix);
	284	for(int iy=1; iy<=ay->GetNbins(); iy++){
	285	rxy[1] = ay->GetBinCenter(iy);
286
287	Double_t rr,ee;
288	pdf.Eval(rxy, rr, ee, kFALSE);
289	if(rr<0. \|\| ee/rr>.15) continue; // 15% relative error
290	//printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
291	h2s->SetBinContent(ix, iy, rr);
292	}
293	}
294
295	// write results to output array
296	//pdf.GetStatus();
297	pdf.Write(Form("%s[%d]", AliPID::ParticleShortName(is), ip));
298	}
299	}
300	delete [] data[0]; delete [] data[1];
301	pdfs->Write();
b91fdd71	302	//fCalib->Close(); delete fCalib;
1ee39b3a	303
	304	return kTRUE; // testing protection
	305	}
	306
	307