[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)
//  - Neural Network (responsible A.Wilk)
//
//   Origin
//   Alex Bercuci  (A.Bercuci@gsi.de)
//
///////////////////////////////////////////////////////////////////////////////

#include <TFile.h>
#include <TROOT.h>
#include <TTree.h>
#include <TMath.h>
#include <TEventList.h>
#include <TH2D.h>
#include <TH2I.h>
#include <TPrincipal.h>
//#include <TVector3.h>
//#include <TLinearFitter.h>
#include <TCanvas.h>
//#include <TEllipse.h>
//#include <TMarker.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

  AliTRDpidRefMaker::CreateOutputObjects();

  // 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));
    for(Int_t is=AliPID::kSPECIES; is--;) {
      h2 = new TH2D(Form("h%s%d", AliPID::ParticleShortName(is), ip), Form("%s ref. dEdx @ 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, 1+ip);
  }
}

/*
//________________________________________________________________________
Float_t* AliTRDpidRefMakerLQ::CookdEdx(AliTRDseedV1 *trklt)
{
// Fill dEdx array for multidim LQ PID

  trklt->CookdEdx(AliTRDpidUtil::kLQslices);
  const Float_t *dedx = trklt->GetdEdx();
  if(dedx[0]+dedx[1] <= 0.) return 0x0;
  if(dedx[2] <= 0.) return 0x0;

  fdEdx[0] = TMath::Log(dedx[0]+dedx[1]);
  fdEdx[1] = TMath::Log(dedx[2]);
  return fdEdx;
}*/

//__________________________________________________________________
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>AliTRDCalPID::kNMom-1){ 
    AliError("Ref fig requested outside definition.");
    return kFALSE;
  }
  if(!gPad){
    AliWarning("Please provide a canvas to draw results.");
    return kFALSE;
  }

  TObjArray *arr = (TObjArray*)fContainer->At(ifig);
  gPad->Divide(3, 2, 1.e-5, 1.e-5); 
  TList *l=gPad->GetListOfPrimitives(); 
  for(Int_t is=0; is<AliPID::kSPECIES; is++){
    ((TVirtualPad*)l->At(is))->cd();
    ((TH2*)arr->At(is))->Draw("cont4z");
  }

  return kTRUE;
}

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


  TFile *fCalib = TFile::Open(Form("TRD.Calib%s.root", GetName()), "update");
  fData = dynamic_cast<TTree*>(gFile->Get(GetName()));
  if (!fData) {
    AliError("Calibration data not available");
    return kFALSE;
  }
  TObjArray *o = 0x0;
  if(!(o = (TObjArray*)gFile->Get(Form("Moni%s", GetName())))){
    AliWarning("Missing monitoring container.");
    return kFALSE;
  }
  fContainer = (TObjArray*)o->Clone("monitor");

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