[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 "../STAT/TKDInterpolator.h"
#include "AliTRDpidRefMakerLQ.h"
#include "Cal/AliTRDCalPID.h"
#include "Cal/AliTRDCalPIDLQ.h"
#include "AliTRDseedV1.h"
#include "AliTRDcalibDB.h"
#include "AliTRDgeometry.h"
#include "info/AliTRDpidInfo.h"

ClassImp(AliTRDpidRefMakerLQ)

//__________________________________________________________________
AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ() 
  : AliTRDpidRefMaker()
  ,fPDF(NULL)
{
  //
  // AliTRDpidRefMakerLQ default constructor
  //
  SetNameTitle("TRDrefMakerLQ", "PID(LQ) Reference Maker");
}

//__________________________________________________________________
AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ(const char *name)
  : AliTRDpidRefMaker(name, "PID(LQ) Reference Maker")
  ,fPDF(NULL)
{
  //
  // AliTRDpidRefMakerLQ default constructor
  //
  DefineOutput(3, TObjArray::Class());
}

//__________________________________________________________________
AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
{
  //
  // AliTRDCalPIDQRef destructor
  //
  if(fPDF){
    //fPDF->Write("PDF_LQ", TObject::kSingleKey);
    fPDF->Delete();
    delete fPDF;
  }
}

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

  fContainer = Histos();
  PostData(1, fContainer);

  //OpenFile(2, "RECREATE");
  fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
  fPDF->SetOwner();fPDF->SetName("PDF_LQ");
  PostData(3, fPDF);
}


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

  if(fContainer) return fContainer;
  fContainer  = new TObjArray(AliTRDCalPID::kNMom);

  // save dE/dx references
  TH1 *h(NULL);
  for(Int_t ip=AliTRDCalPID::kNMom; ip--;){ 
    TObjArray *arr = new TObjArray(2*AliPID::kSPECIES);
    arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
    for(Int_t is=AliPID::kSPECIES; is--;) {
      h = new TH1D(Form("h1%s%d", AliPID::ParticleShortName(is), ip), Form("1D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12.);
      h->GetXaxis()->SetTitle("log(dE/dx) [au]");
      h->GetYaxis()->SetTitle("#");
      h->SetLineColor(AliTRDCalPIDLQ::GetPartColor(is));
      arr->AddAt(h, is);
    }
    for(Int_t is=AliPID::kSPECIES; is--;) {
      h = new TH2D(Form("h2%s%d", AliPID::ParticleShortName(is), ip), Form("2D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12., 50, 6.5, 11.);
      h->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
      h->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
      h->GetZaxis()->SetTitle("#");
      arr->AddAt(h, AliPID::kSPECIES+is);
    }
    fContainer->AddAt(arr, ip);
  }
  fNRefFigures=AliTRDCalPID::kNMom;
  return fContainer;
}


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

  if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
    AliError(Form("File TRD.Calib%s.root not readable", GetName()));
    return NULL;
  }
  AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
  cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
  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);TH1 *h(NULL);
  if(!(arr = (TObjArray*)fContainer->At(ifig))){
    AliError(Form("PDF container @ pBin[%d] missing.", ifig));
    return kFALSE;
  }
  gPad->Divide(5, 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(!(h=(TH1*)arr->At(is))){
      AliError(Form("1D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
      return kFALSE;
    }
    h->GetYaxis()->SetRangeUser(0., 1.2);
    h->Draw("l");

    ((TVirtualPad*)l->At(AliPID::kSPECIES+is))->cd();
    if(!(h=(TH1*)arr->At(AliPID::kSPECIES+is))){
      AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
      return kFALSE;
    }
    h->Draw("cont4z");
  }

  return kTRUE;
}


//________________________________________________________________________
Bool_t AliTRDpidRefMakerLQ::Load(const Char_t *file, const Char_t *dir)
{
// Load tree with data in case of detached PostProcess processing. 

  if(gSystem->AccessPathName(file, kReadPermission)){
    AliError(Form("File %s not readable", file));
    return kFALSE;
  }
  if(!TFile::Open(file)) {
    AliError(Form("File %s corrupted", file));
    return kFALSE;
  }
  if(!gFile->cd(dir)){
    AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
    return kFALSE;
  }
  if (!(fData = dynamic_cast<TTree*>(gDirectory->Get("PIDrefMaker")))) {
    AliError("PIDref Tree not available");
    return kFALSE;
  }
  LinkPIDdata();

  TObjArray *o(NULL);
/*  if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
    AliWarning(Form("Monitor container Moni%s not available.", name));
    return kFALSE;
  }
  fContainer = (TObjArray*)o->Clone("monitor");
  fNRefFigures=AliTRDCalPID::kNMom;
*/
  // temporary until new calibration data are being produced
  Histos();


  if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
    AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
    return kFALSE;
  }
  if(!(o = (TObjArray*)gFile->Get(Form("PDF_LQ")))) {
    AliInfo("PDF container not available. Create.");
    fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
    fPDF->SetOwner();fPDF->SetName("PDF_LQ");
  } else fPDF = (TObjArray*)o->Clone("PDF_LQ");

  return kTRUE;
}


//________________________________________________________________________
void AliTRDpidRefMakerLQ::UserExec(Option_t */*opt*/)
{
// Load PID data into local data storage

  if(!(fInfo   = dynamic_cast<TObjArray*>(GetInputData(3)))) return;

  AliDebug(2, Form("ENTRIES pid[%d]\n", fInfo->GetEntriesFast()));
  AliTRDpidInfo *pid(NULL);
  const AliTRDpidInfo::AliTRDpidData *data(NULL);
  Char_t s(-1);
  for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
    if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
    if((s=pid->GetPID())<0) continue;
    for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
      data=pid->GetData(itrklt);
      Int_t ip(data->Momentum());

      Double_t dedx[] = {0., 0.};
      if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
      ((TH2*)((TObjArray*)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]);
      AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE);
      ((TH1*)((TObjArray*)fContainer->At(ip))->At(s))->Fill(dedx[0]+dedx[1]);
    }
  }
}


//________________________________________________________________________
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
// 

  TCanvas *fMonitor(NULL);
  // allocate working storage
  const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
  Float_t *data[2*kWS], *data1D[kWS];
  for(Int_t i(0); i<kWS; i++){ 
    data1D[i]   = new Float_t[kMaxStat];
    data[i]     = new Float_t[kMaxStat];
    data[kWS+i] = new Float_t[kMaxStat];
  }
  Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));

  AliDebug(1, Form("Loading %d entries.", (Int_t)fData->GetEntries()));
  for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
    if(!(fData->GetEntry(itrk))) continue;
    Int_t sbin(fPIDdataArray->GetPID());
    for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
      Int_t pbin(fPIDdataArray->GetData(itrklt)->Momentum());
      
      Double_t dedx[] = {0., 0.}, 
               dedx1D[] = {0., 0.};
      if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx)) continue;
      AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx1D, kFALSE);
      Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
      if(ndata[idx]==kMaxStat) continue;

      // store data
      data1D[idx][ndata[idx]]   = dedx1D[0];
      data[idx][ndata[idx]]     = dedx[0];
      data[idx+kWS][ndata[idx]] = dedx[1];
      ndata[idx]++;
    }
  }

  TKDPDF *pdf(NULL);
  Int_t in(0); Float_t par[6], *pp(NULL);
  for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){ 
    for(Int_t is=AliPID::kSPECIES; is--;) {
      // estimate bucket statistics
      Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
            nb(kMinBuckets), // number of buckets
            ns((Int_t)(((Float_t)(ndata[idx]))/nb));    //statistics/bucket
            
      AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
      if(ns<Int_t(kMinStat)){
        AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
        continue;
      }

      // build helper 1D PDF
      pdf = new TKDPDF(ndata[idx], 1, ns, &data1D[idx]);
      pdf->SetCOG(kFALSE);
      pdf->SetWeights();
      //pdf->SetStore();
      pdf->SetAlpha(15.);
      //pdf.GetStatus();
      fPDF->AddAt(pdf, idx);
      in=pdf->GetNTNodes(); pp=&par[0];
      while(in--){
        const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
        nn->GetCOG(pp);
        Double_t p(par[0]),r,e;
        pdf->Eval(&p,r,e,1);
      }

      // build helper 2D PDF
      Float_t *ldata[2]={data[idx], data[kWS+idx]};
      pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
      pdf->SetCOG(kFALSE);
      pdf->SetWeights();
      //pdf->SetStore();
      pdf->SetAlpha(5.);
      //pdf.GetStatus();
      fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2));
      in=pdf->GetNTNodes(); pp=&par[0];
      while(in--){
        const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
        nn->GetCOG(pp);
        Double_t p[] = {par[0], par[1]}, r,e;
        pdf->Eval(p,r,e,1);
      }
/*
      Int_t nnodes = pdf.GetNTNodes(),
            nside = Int_t(0.05*nnodes),
            nzeros = 4*(nside+1);
      printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
    

      // Build interpolator on the pdf skeleton
      TKDInterpolator interpolator(2, nnodes+nzeros); 
      for(Int_t in=nnodes; in--;)
        interpolator.SetNode(in, *pdf.GetNodeInfo(in));
      TKDNodeInfo *nodes = new TKDNodeInfo[nzeros], *node = &nodes[0];
      Float_t ax0min, ax0max, ax1min, ax1max;
      pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
      pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
      printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);

      Int_t jn = nnodes; 
      SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
      SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
      SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
      SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
      delete [] nodes;
      Int_t in=nnodes; Float_t par[6], *pp=&par[0];
      while(in--){
        const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
        nn->GetCOG(pp);
        //printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
        Double_t p[] = {par[0], par[1]}, r,e;
        interpolator.Eval(p,r,e,1);
      }
*/

      // visual on-line monitoring
      if(HasOnlineMonitor()){
        if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
        pdf->DrawProjection();
        fMonitor->Modified(); fMonitor->Update(); 
        fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
      }

      //fContainer->ls();
      // save a discretization of the PDF for result monitoring
      Double_t rxy[]={0.,0.};
      TH2 *h2s = (TH2D*)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+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);
        }
      }

      if(!(pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx)))){
        AliWarning(Form("Missing pdf for model id[%d]", idx));
        continue;
      }
      TH1 *h1 = (TH1D*)((TObjArray*)fContainer->At(ip))->At(is);
      ax = h1->GetXaxis();
      h1->Clear();
      for(int ix=1; ix<=ax->GetNbins(); ix++){
        rxy[0] = ax->GetBinCenter(ix);
      
        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);
        h1->SetBinContent(ix, rr);
      }
    }
  }
  for(Int_t i(0); i<kWS; i++){ 
    delete [] data1D[i];
    delete [] data[i];
    delete [] data[i+kWS];
  }
  return kTRUE;
}


//__________________________________________________________________
void AliTRDpidRefMakerLQ::SetZeroes(TKDInterpolator *interpolator, TKDNodeInfo *node, Int_t n, Int_t& idx, Float_t x, Float_t dx, Float_t y, Float_t dy, const Char_t opt)
{
// Set extra nodes to ensure boundary conditions
  
  printf("SetZeroes(%c)\n", opt);
  Float_t par[6], val[] = {0., 1.};
  Int_t a[6];
  if(opt=='x'){
    a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
  } else if(opt=='y'){
    a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
  } else return;
  Float_t tmp;
  par[a[1]] = y;
  par[a[4]] = y; par[a[5]] = y+dy;
  if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
  for(Int_t in=n; in--; node++, idx++, x+=dx){
    par[a[0]] = x+.5*dx;
    par[a[2]] = x;  par[a[3]] = x+dx;
    if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
    node->SetNode(2, par, val);
    printf("\n\tnode[%d]\n", idx); node->Print();
    interpolator->SetNode(idx, *node);
  }
  par[a[0]] = x;
  par[a[2]] = x;  par[a[3]] = x+dx;
  if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
  node->SetNode(2, par, val);
  printf("\n\tnode[%d]\n", idx); node->Print();
  interpolator->SetNode(idx, *node);node++;idx++;
}
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"
4826d5b1	42	#include "../STAT/TKDInterpolator.h"
1ee39b3a	43	#include "AliTRDpidRefMakerLQ.h"
da27d983	44	#include "Cal/AliTRDCalPID.h"
da27d983	45	#include "Cal/AliTRDCalPIDLQ.h"
1ee39b3a	46	#include "AliTRDseedV1.h"
	47	#include "AliTRDcalibDB.h"
	48	#include "AliTRDgeometry.h"
4826d5b1	49	#include "info/AliTRDpidInfo.h"
1ee39b3a	50
	51	ClassImp(AliTRDpidRefMakerLQ)
	52
	53	//__________________________________________________________________
4826d5b1	54	AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ()
705f8b0a	55	: AliTRDpidRefMaker()
a5a3321d	56	,fPDF(NULL)
1ee39b3a	57	{
	58	//
	59	// AliTRDpidRefMakerLQ default constructor
	60	//
4fa7d600	61	SetNameTitle("TRDrefMakerLQ", "PID(LQ) Reference Maker");
705f8b0a	62	}
	63
	64	//__________________________________________________________________
	65	AliTRDpidRefMakerLQ::AliTRDpidRefMakerLQ(const char *name)
	66	: AliTRDpidRefMaker(name, "PID(LQ) Reference Maker")
a5a3321d	67	,fPDF(NULL)
705f8b0a	68	{
	69	//
	70	// AliTRDpidRefMakerLQ default constructor
	71	//
	72	DefineOutput(3, TObjArray::Class());
1ee39b3a	73	}
	74
	75	//__________________________________________________________________
	76	AliTRDpidRefMakerLQ::~AliTRDpidRefMakerLQ()
	77	{
	78	//
	79	// AliTRDCalPIDQRef destructor
	80	//
4826d5b1	81	if(fPDF){
705f8b0a	82	//fPDF->Write("PDF_LQ", TObject::kSingleKey);
4826d5b1	83	fPDF->Delete();
	84	delete fPDF;
	85	}
1ee39b3a	86	}
1ee39b3a	87
b91fdd71	88	// //________________________________________________________________________
f8f46e4d	89	void AliTRDpidRefMakerLQ::UserCreateOutputObjects()
1ee39b3a	90	{
	91	// Create histograms
	92	// Called once
	93
b91fdd71	94	fContainer = Histos();
068e2c00	95	PostData(1, fContainer);
4826d5b1	96
a5a3321d	97	//OpenFile(2, "RECREATE");
092fbd8a	98	fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
092fbd8a	99	fPDF->SetOwner();fPDF->SetName("PDF_LQ");
068e2c00	100	PostData(3, fPDF);
b91fdd71	101	}
	102
	103
	104	//__________________________________________________________________
	105	TObjArray* AliTRDpidRefMakerLQ::Histos()
	106	{
	107	// Create histograms
	108
a5a3321d	109	if(fContainer) return fContainer;
a5a3321d	110	fContainer = new TObjArray(AliTRDCalPID::kNMom);
1ee39b3a	111
1ee39b3a	112	// save dE/dx references
092fbd8a	113	TH1 *h(NULL);
1ee39b3a	114	for(Int_t ip=AliTRDCalPID::kNMom; ip--;){
092fbd8a	115	TObjArray arr = new TObjArray(2AliPID::kSPECIES);
b91fdd71	116	arr->SetName(Form("Pbin%02d", ip)); arr->SetOwner();
1ee39b3a	117	for(Int_t is=AliPID::kSPECIES; is--;) {
092fbd8a	118	h = new TH1D(Form("h1%s%d", AliPID::ParticleShortName(is), ip), Form("1D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12.);
	119	h->GetXaxis()->SetTitle("log(dE/dx) [au]");
	120	h->GetYaxis()->SetTitle("#");
	121	h->SetLineColor(AliTRDCalPIDLQ::GetPartColor(is));
	122	arr->AddAt(h, is);
	123	}
	124	for(Int_t is=AliPID::kSPECIES; is--;) {
	125	h = new TH2D(Form("h2%s%d", AliPID::ParticleShortName(is), ip), Form("2D %s @ Pbin[%d]", AliPID::ParticleName(is), ip), 50, 7., 12., 50, 6.5, 11.);
	126	h->GetXaxis()->SetTitle("log(dE/dx_{am}) [au]");
	127	h->GetYaxis()->SetTitle("log(dE/dx_{dr}) [au]");
	128	h->GetZaxis()->SetTitle("#");
	129	arr->AddAt(h, AliPID::kSPECIES+is);
1ee39b3a	130	}
a5a3321d	131	fContainer->AddAt(arr, ip);
1ee39b3a	132	}
a5a3321d	133	fNRefFigures=AliTRDCalPID::kNMom;
b91fdd71	134	return fContainer;
1ee39b3a	135	}
1ee39b3a	136
1ee39b3a	137
1ee39b3a	138	//__________________________________________________________________
4826d5b1	139	TObject* AliTRDpidRefMakerLQ::GetOCDBEntry(Option_t /opt*/)
1ee39b3a	140	{
	141	// Steer loading of OCDB LQ PID
	142
4826d5b1	143	if(gSystem->AccessPathName(Form("TRD.Calib%s.root", GetName()), kReadPermission)){
	144	AliError(Form("File TRD.Calib%s.root not readable", GetName()));
	145	return NULL;
	146	}
1ee39b3a	147	AliTRDCalPIDLQ *cal = new AliTRDCalPIDLQ("pidLQ", "LQ TRD PID object");
4826d5b1	148	cal->LoadReferences(Form("TRD.Calib%s.root", GetName()));
1ee39b3a	149	return cal;
	150	}
	151
	152	//__________________________________________________________________
	153	Bool_t AliTRDpidRefMakerLQ::GetRefFigure(Int_t ifig)
	154	{
	155	// Steering reference picture
	156
b91fdd71	157	if(ifig<0 \|\| ifig>=fNRefFigures){
1ee39b3a	158	AliError("Ref fig requested outside definition.");
	159	return kFALSE;
	160	}
	161	if(!gPad){
	162	AliWarning("Please provide a canvas to draw results.");
	163	return kFALSE;
	164	}
	165
092fbd8a	166	TObjArray arr(NULL);TList l(NULL);TH1 *h(NULL);
	167	if(!(arr = (TObjArray*)fContainer->At(ifig))){
	168	AliError(Form("PDF container @ pBin[%d] missing.", ifig));
	169	return kFALSE;
	170	}
	171	gPad->Divide(5, 2, 1.e-5, 1.e-5);l=gPad->GetListOfPrimitives();
	172	for(Int_t is=0; is<AliPID::kSPECIES; is++){
	173	((TVirtualPad*)l->At(is))->cd();
	174	if(!(h=(TH1*)arr->At(is))){
	175	AliError(Form("1D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
b91fdd71	176	return kFALSE;
b91fdd71	177	}
092fbd8a	178	h->GetYaxis()->SetRangeUser(0., 1.2);
	179	h->Draw("l");
	180
	181	((TVirtualPad*)l->At(AliPID::kSPECIES+is))->cd();
	182	if(!(h=(TH1*)arr->At(AliPID::kSPECIES+is))){
	183	AliError(Form("2D for %s @ pBin[%d] missing.", AliPID::ParticleShortName(is), ifig));
b91fdd71	184	return kFALSE;
b91fdd71	185	}
092fbd8a	186	h->Draw("cont4z");
b91fdd71	187	}
092fbd8a	188
b91fdd71	189	return kTRUE;
	190	}
	191
	192
	193	//________________________________________________________________________
fe1d1beb	194	Bool_t AliTRDpidRefMakerLQ::Load(const Char_t file, const Char_t dir)
b91fdd71	195	{
64d57299	196	// Load tree with data in case of detached PostProcess processing.
64d57299	197
fe1d1beb	198	if(gSystem->AccessPathName(file, kReadPermission)){
fe1d1beb	199	AliError(Form("File %s not readable", file));
b91fdd71	200	return kFALSE;
b91fdd71	201	}
fe1d1beb	202	if(!TFile::Open(file)) {
fe1d1beb	203	AliError(Form("File %s corrupted", file));
b91fdd71	204	return kFALSE;
b91fdd71	205	}
fe1d1beb	206	if(!gFile->cd(dir)){
	207	AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
	208	return kFALSE;
	209	}
	210	if (!(fData = dynamic_cast<TTree*>(gDirectory->Get("PIDrefMaker")))) {
	211	AliError("PIDref Tree not available");
b91fdd71	212	return kFALSE;
1ee39b3a	213	}
b91fdd71	214	LinkPIDdata();
1ee39b3a	215
b91fdd71	216	TObjArray *o(NULL);
092fbd8a	217	/* if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
b91fdd71	218	AliWarning(Form("Monitor container Moni%s not available.", name));
	219	return kFALSE;
	220	}
	221	fContainer = (TObjArray*)o->Clone("monitor");
4826d5b1	222	fNRefFigures=AliTRDCalPID::kNMom;
092fbd8a	223	*/
	224	// temporary until new calibration data are being produced
	225	Histos();
4826d5b1	226
	227
	228	if(!TFile::Open(Form("TRD.Calib%s.root", GetName()), "UPDATE")){
	229	AliError(Form("File TRD.Calib%s.root corrupted", GetName()));
	230	return kFALSE;
	231	}
092fbd8a	232	if(!(o = (TObjArray*)gFile->Get(Form("PDF_LQ")))) {
	233	AliInfo("PDF container not available. Create.");
	234	fPDF = new TObjArray(AliTRDCalPIDLQ::GetNrefs());
	235	fPDF->SetOwner();fPDF->SetName("PDF_LQ");
	236	} else fPDF = (TObjArray*)o->Clone("PDF_LQ");
4826d5b1	237
1ee39b3a	238	return kTRUE;
	239	}
	240
4826d5b1	241
b91fdd71	242	//________________________________________________________________________
f8f46e4d	243	void AliTRDpidRefMakerLQ::UserExec(Option_t /opt*/)
b91fdd71	244	{
4826d5b1	245	// Load PID data into local data storage
4826d5b1	246
87c9fc00	247	if(!(fInfo = dynamic_cast<TObjArray*>(GetInputData(3)))) return;
87c9fc00	248
a5a3321d	249	AliDebug(2, Form("ENTRIES pid[%d]\n", fInfo->GetEntriesFast()));
4826d5b1	250	AliTRDpidInfo *pid(NULL);
	251	const AliTRDpidInfo::AliTRDpidData *data(NULL);
	252	Char_t s(-1);
	253	for(Int_t itrk=fInfo->GetEntriesFast(); itrk--;){
	254	if(!(pid=(AliTRDpidInfo*)fInfo->At(itrk))) continue;
	255	if((s=pid->GetPID())<0) continue;
	256	for(Int_t itrklt=pid->GetNtracklets();itrklt--;){
	257	data=pid->GetData(itrklt);
	258	Int_t ip(data->Momentum());
4826d5b1	259
	260	Double_t dedx[] = {0., 0.};
	261	if(!AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx)) continue;
092fbd8a	262	((TH2)((TObjArray)fContainer->At(ip))->At(s+Int_t(AliPID::kSPECIES)))->Fill(dedx[0], dedx[1]);
092fbd8a	263	AliTRDCalPIDLQ::CookdEdx(data->fdEdx, dedx, kFALSE);
a5a3321d	264	((TH1)((TObjArray)fContainer->At(ip))->At(s))->Fill(dedx[0]+dedx[1]);
4826d5b1	265	}
4826d5b1	266	}
b91fdd71	267	}
	268
	269
1ee39b3a	270	//________________________________________________________________________
	271	Bool_t AliTRDpidRefMakerLQ::PostProcess()
	272	{
	273	// Analyse merged dedx = f(p) distributions.
	274	// - select momentum - species bins
	275	// - rotate to principal components
	276	// - locally interpolate with TKDPDF
	277	// - save interpolation to monitoring histograms
	278	// - write pdf to file for loading to OCDB
	279	//
	280
674c7d26	281	TCanvas *fMonitor(NULL);
b91fdd71	282	// allocate working storage
674c7d26	283	const Int_t kWS(AliPID::kSPECIES*AliTRDCalPID::kNMom);
092fbd8a	284	Float_t data[2kWS], *data1D[kWS];
	285	for(Int_t i(0); i<kWS; i++){
	286	data1D[i] = new Float_t[kMaxStat];
	287	data[i] = new Float_t[kMaxStat];
	288	data[kWS+i] = new Float_t[kMaxStat];
	289	}
674c7d26	290	Int_t ndata[kWS]; memset(ndata, 0, kWS*sizeof(Int_t));
674c7d26	291
ca50a37e	292	AliDebug(1, Form("Loading %d entries.", (Int_t)fData->GetEntries()));
674c7d26	293	for(Int_t itrk=0; itrk < fData->GetEntries(); itrk++){
674c7d26	294	if(!(fData->GetEntry(itrk))) continue;
a5a3321d	295	Int_t sbin(fPIDdataArray->GetPID());
	296	for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
	297	Int_t pbin(fPIDdataArray->GetData(itrklt)->Momentum());
674c7d26	298
092fbd8a	299	Double_t dedx[] = {0., 0.},
092fbd8a	300	dedx1D[] = {0., 0.};
a5a3321d	301	if(!AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx)) continue;
a5a3321d	302	AliTRDCalPIDLQ::CookdEdx(fPIDdataArray->GetData(itrklt)->fdEdx, dedx1D, kFALSE);
674c7d26	303	Int_t idx=AliTRDCalPIDLQ::GetModelID(pbin,sbin);
	304	if(ndata[idx]==kMaxStat) continue;
	305
	306	// store data
092fbd8a	307	data1D[idx][ndata[idx]] = dedx1D[0];
092fbd8a	308	data[idx][ndata[idx]] = dedx[0];
674c7d26	309	data[idx+kWS][ndata[idx]] = dedx[1];
	310	ndata[idx]++;
	311	}
	312	}
092fbd8a	313
	314	TKDPDF *pdf(NULL);
	315	Int_t in(0); Float_t par[6], *pp(NULL);
4826d5b1	316	for(Int_t ip=0;ip<AliTRDCalPID::kNMom;ip++){
1ee39b3a	317	for(Int_t is=AliPID::kSPECIES; is--;) {
1ee39b3a	318	// estimate bucket statistics
674c7d26	319	Int_t idx(AliTRDCalPIDLQ::GetModelID(ip,is)),
674c7d26	320	nb(kMinBuckets), // number of buckets
7fe4e88b	321	ns((Int_t)(((Float_t)(ndata[idx]))/nb)); //statistics/bucket
1ee39b3a	322
674c7d26	323	AliDebug(2, Form("pBin[%d] sBin[%d] n[%d] ns[%d] nb[%d]", ip, is, ndata[idx], ns, nb));
1ee39b3a	324	if(ns<Int_t(kMinStat)){
674c7d26	325	AliWarning(Form("Not enough entries [%d] for %s[%d].", ndata[idx], AliPID::ParticleShortName(is), ip));
1ee39b3a	326	continue;
	327	}
	328
092fbd8a	329	// build helper 1D PDF
	330	pdf = new TKDPDF(ndata[idx], 1, ns, &data1D[idx]);
	331	pdf->SetCOG(kFALSE);
	332	pdf->SetWeights();
	333	//pdf->SetStore();
	334	pdf->SetAlpha(15.);
	335	//pdf.GetStatus();
	336	fPDF->AddAt(pdf, idx);
	337	in=pdf->GetNTNodes(); pp=&par[0];
	338	while(in--){
	339	const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
	340	nn->GetCOG(pp);
	341	Double_t p(par[0]),r,e;
	342	pdf->Eval(&p,r,e,1);
	343	}
	344
	345	// build helper 2D PDF
674c7d26	346	Float_t *ldata[2]={data[idx], data[kWS+idx]};
092fbd8a	347	pdf = new TKDPDF(ndata[idx], 2, ns, ldata);
4826d5b1	348	pdf->SetCOG(kFALSE);
4826d5b1	349	pdf->SetWeights();
674c7d26	350	//pdf->SetStore();
4826d5b1	351	pdf->SetAlpha(5.);
4826d5b1	352	//pdf.GetStatus();
092fbd8a	353	fPDF->AddAt(pdf, AliTRDCalPIDLQ::GetModelID(ip,is, 2));
092fbd8a	354	in=pdf->GetNTNodes(); pp=&par[0];
4826d5b1	355	while(in--){
	356	const TKDNodeInfo *nn = pdf->GetNodeInfo(in);
	357	nn->GetCOG(pp);
4826d5b1	358	Double_t p[] = {par[0], par[1]}, r,e;
4826d5b1	359	pdf->Eval(p,r,e,1);
1ee39b3a	360	}
4826d5b1	361	/*
	362	Int_t nnodes = pdf.GetNTNodes(),
	363	nside = Int_t(0.05*nnodes),
	364	nzeros = 4*(nside+1);
	365	printf("nnodes[%d] nside[%d] nzeros[%d]\n", nnodes, nside, nzeros);
	366
	367
	368	// Build interpolator on the pdf skeleton
	369	TKDInterpolator interpolator(2, nnodes+nzeros);
	370	for(Int_t in=nnodes; in--;)
	371	interpolator.SetNode(in, *pdf.GetNodeInfo(in));
	372	TKDNodeInfo nodes = new TKDNodeInfo[nzeros], node = &nodes[0];
	373	Float_t ax0min, ax0max, ax1min, ax1max;
	374	pdf.GetRange(0, ax0min, ax0max); Float_t dx = (ax0max-ax0min)/nside;
	375	pdf.GetRange(1, ax1min, ax1max); Float_t dy = (ax1max-ax1min)/nside;
	376	printf("x=[%f %f] y[%f %f]\n", ax0min, ax0max, ax1min, ax1max);
	377
	378	Int_t jn = nnodes;
	379	SetZeroes(&interpolator, node, nside, jn, ax0min, dx, ax1min, -dy, 'x');
	380	SetZeroes(&interpolator, node, nside, jn, ax1min, dy, ax0max, dx, 'y');
	381	SetZeroes(&interpolator, node, nside, jn, ax0max,-dx, ax1max, dy, 'x');
	382	SetZeroes(&interpolator, node, nside, jn ,ax1max, -dy, ax0min, -dx, 'y');
	383	delete [] nodes;
674c7d26	384	Int_t in=nnodes; Float_t par[6], *pp=&par[0];
4826d5b1	385	while(in--){
	386	const TKDNodeInfo *nn = interpolator.GetNodeInfo(in);
	387	nn->GetCOG(pp);
	388	//printf("evaluate for node[%d] @ [%f %f]\n",in, par[0], par[1]);
	389	Double_t p[] = {par[0], par[1]}, r,e;
	390	interpolator.Eval(p,r,e,1);
	391	}
	392	*/
092fbd8a	393
4826d5b1	394	// visual on-line monitoring
674c7d26	395	if(HasOnlineMonitor()){
	396	if(!fMonitor) fMonitor = new TCanvas("cc", "PDF 2D LQ", 500, 500);
	397	pdf->DrawProjection();
	398	fMonitor->Modified(); fMonitor->Update();
	399	fMonitor->SaveAs(Form("pdf_%s%02d.gif", AliPID::ParticleShortName(is), ip));
	400	}
1ee39b3a	401
4826d5b1	402	//fContainer->ls();
1ee39b3a	403	// save a discretization of the PDF for result monitoring
4826d5b1	404	Double_t rxy[]={0.,0.};
092fbd8a	405	TH2 h2s = (TH2D)((TObjArray*)fContainer->At(ip))->At(AliPID::kSPECIES+is);
1ee39b3a	406	TAxis ax = h2s->GetXaxis(), ay = h2s->GetYaxis();
	407	h2s->Clear();
	408	for(int ix=1; ix<=ax->GetNbins(); ix++){
	409	rxy[0] = ax->GetBinCenter(ix);
	410	for(int iy=1; iy<=ay->GetNbins(); iy++){
	411	rxy[1] = ay->GetBinCenter(iy);
	412
	413	Double_t rr,ee;
4826d5b1	414	pdf->Eval(rxy, rr, ee, kFALSE);
1ee39b3a	415	if(rr<0. \|\| ee/rr>.15) continue; // 15% relative error
	416	//printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
	417	h2s->SetBinContent(ix, iy, rr);
	418	}
	419	}
092fbd8a	420
7fe4e88b	421	if(!(pdf=dynamic_cast<TKDPDF*>(fPDF->At(idx)))){
	422	AliWarning(Form("Missing pdf for model id[%d]", idx));
	423	continue;
	424	}
092fbd8a	425	TH1 h1 = (TH1D)((TObjArray*)fContainer->At(ip))->At(is);
	426	ax = h1->GetXaxis();
	427	h1->Clear();
	428	for(int ix=1; ix<=ax->GetNbins(); ix++){
	429	rxy[0] = ax->GetBinCenter(ix);
	430
	431	Double_t rr,ee;
	432	pdf->Eval(rxy, rr, ee, kFALSE);
	433	if(rr<0. \|\| ee/rr>.15) continue; // 15% relative error
	434	//printf("x[%2d] x[%2d] r[%f] e[%f]\n", ix, iy, rr, ee);
	435	h1->SetBinContent(ix, rr);
	436	}
1ee39b3a	437	}
1ee39b3a	438	}
092fbd8a	439	for(Int_t i(0); i<kWS; i++){
	440	delete [] data1D[i];
	441	delete [] data[i];
	442	delete [] data[i+kWS];
	443	}
4826d5b1	444	return kTRUE;
1ee39b3a	445	}
	446
	447
4826d5b1	448	//__________________________________________________________________
	449	void AliTRDpidRefMakerLQ::SetZeroes(TKDInterpolator interpolator, TKDNodeInfo node, Int_t n, Int_t& idx, Float_t x, Float_t dx, Float_t y, Float_t dy, const Char_t opt)
	450	{
	451	// Set extra nodes to ensure boundary conditions
	452
	453	printf("SetZeroes(%c)\n", opt);
	454	Float_t par[6], val[] = {0., 1.};
	455	Int_t a[6];
	456	if(opt=='x'){
	457	a[0]=0; a[1]=1; a[2]=2; a[3]=3; a[4]=4; a[5]=5;
	458	} else if(opt=='y'){
	459	a[0]=1; a[1]=0; a[2]=4; a[3]=5; a[4]=2; a[5]=3;
	460	} else return;
	461	Float_t tmp;
	462	par[a[1]] = y;
	463	par[a[4]] = y; par[a[5]] = y+dy;
	464	if(dy<0.){tmp=par[a[4]]; par[a[4]]=par[a[5]]; par[a[5]]=tmp;}
	465	for(Int_t in=n; in--; node++, idx++, x+=dx){
	466	par[a[0]] = x+.5*dx;
	467	par[a[2]] = x; par[a[3]] = x+dx;
	468	if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
	469	node->SetNode(2, par, val);
	470	printf("\n\tnode[%d]\n", idx); node->Print();
	471	interpolator->SetNode(idx, *node);
	472	}
	473	par[a[0]] = x;
	474	par[a[2]] = x; par[a[3]] = x+dx;
	475	if(dx<0.){tmp=par[a[2]]; par[a[2]]=par[a[3]]; par[a[3]]=tmp;}
	476	node->SetNode(2, par, val);
	477	printf("\n\tnode[%d]\n", idx); node->Print();
	478	interpolator->SetNode(idx, *node);node++;idx++;
	479	}
	480