[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPIDLQ.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$ */

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Container for the distributions of dE/dx and the time bin of the     //
// max. cluster for electrons and pions                                 //
//                                                                      //
// Authors:                                                             //
//   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
//   Alex Bercuci <a.bercuci@gsi.de>                                    //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include <TH2F.h>
#include <TFile.h>
#include <TROOT.h>

#include "AliLog.h"
#include "AliPID.h"

#include "../../STAT/TKDPDF.h"
#include "AliTRDCalPIDLQ.h"
#include "AliTRDcalibDB.h"

ClassImp(AliTRDCalPIDLQ)

Float_t AliTRDCalPIDLQ::fgTrackSegLength[kNLength] = { 3.7, 3.9, 4.2, 5.0 };

//_________________________________________________________________________
AliTRDCalPIDLQ::AliTRDCalPIDLQ()
  :AliTRDCalPID("pid", "LQ PID references for TRD")
{
  //
  //  The Default constructor
  //

  //Init();

}

//_________________________________________________________________________
AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
  :AliTRDCalPID(name,title)
{
  //
  //  The main constructor
  //
  
  Init();

}

//_________________________________________________________________________
AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
{
  //
  // Destructor
  //
  
}

//_________________________________________________________________________
Bool_t AliTRDCalPIDLQ::LoadPDF(TDirectoryFile *d)
{
  // Read histograms
  TKDPDF *pdf = 0x0;
  for (Int_t is=0; is < AliPID::kSPECIES; is++){
    for (Int_t ip = 0; ip < kNMom; ip++){
      if(!(pdf = (TKDPDF*)d->Get(Form("%s[%d]", AliPID::ParticleShortName(is), ip)))){ 
        AliWarning(Form("Reference for %s[%d] missing.", AliPID::ParticleShortName(is), ip));
        continue;
      }
      fModel->AddAt(pdf->Clone(), GetModelID(ip, is, 0));
    }
  }
  return kTRUE;
}
  //
//_________________________________________________________________________
Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
{
  //
  // Read the TRD dEdx histograms.
  //

  Int_t nTimeBins = 22;
  // Get number of time bins from CDB
  AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
  if(!calibration){
    AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
  }else{
    if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
    else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
  }

  
  // Read histogram Root file  
  TFile *histFile = TFile::Open(refFile, "READ");
  if (!histFile || !histFile->IsOpen()) {
    AliError(Form("Opening TRD histgram file %s failed", refFile));
    return kFALSE;
  }
  gROOT->cd();

  // Read histograms
  for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
    for (Int_t imom = 0; imom < kNMom; imom++){
      TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/*, ilength*/))->Clone();
      hist->Scale(1./hist->Integral());
      fModel->AddAt(hist, GetModelID(imom, iparticle, 0));

// 			if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
// 
// 			TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fPartSymb[iparticle], imom))->Clone();
// 			if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fPartSymb[iparticle], imom, nTimeBins));
// 			ht->Scale(1./ht->Integral());
// 			fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
    }
  }
  
  histFile->Close();
  delete histFile;
  
  // Number of bins and bin size
  //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
  //fNbins   = hist->GetNbinsX();
  //fBinSize = hist->GetBinWidth(1);
  
  return kTRUE;


}

//_________________________________________________________________________
TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const          // iType not needed
{
  //
  // Returns one selected dEdx histogram
  //

  if (iType < 0 || iType >= AliPID::kSPECIES) return 0x0;
  if(ip<0 || ip>= kNMom ) return 0x0;

  AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fgTrackMomentum[ip]));
  
  return fModel->At(GetModelID(ip, iType, iplane));
}

//_________________________________________________________________________
Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
                                      , const Float_t * const dedx
                                      , Float_t length, Int_t iplane) const
{
  //
  // Core function of AliTRDCalPID class for calculating the
  // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
  // given momentum "mom" and a given dE/dx (slice "dedx") yield per
  // layer
  //

  if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
    
  //Double_t dedx   = dedx1/fMeanChargeRatio;
  
  // find the interval in momentum and track segment length which applies for this data
  Int_t ilength = 1;
  while(ilength<kNLength-1 && length>fgTrackSegLength[ilength]){
    ilength++;
  }
  Int_t imom = 1;
  while(imom<kNMom-1 && mom>fgTrackMomentum[imom]) imom++;

  Int_t nbinsx, nbinsy;
  TAxis *ax = 0x0, *ay = 0x0;
  Double_t lq1, lq2;
  Double_t mom1 = fgTrackMomentum[imom-1], mom2 = fgTrackMomentum[imom];
  TH2 *hist = 0x0;
  if(!(hist = (TH2D*)fModel->At(GetModelID(imom-1, spec, iplane)))){
    AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
    AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
    return 0.;
  }
  ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
  ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
  Float_t x = dedx[0]+dedx[1], y = dedx[2];
        Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
  Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
  if(kX)
    if(kY) lq1 = hist->GetBinContent( hist->FindBin(x, y)); 
    //fEstimator->Estimate2D2(hist, x, y);
    else lq1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
  else
    if(kY) lq1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
    else lq1 = hist->GetBinContent(nbinsx, nbinsy);


  if(!(hist = (TH2D*)fModel->At(GetModelID(imom, spec, iplane)))){
    AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
    AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom, spec, iplane)));
    return lq1;
  }
  if(kX)
    if(kY) lq2 = hist->GetBinContent( hist->FindBin(x, y)); 
    //fEstimator->Estimate2D2(hist, x, y);
    else lq2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
  else
    if(kY) lq2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
    else lq2 = hist->GetBinContent(nbinsx, nbinsy);
  
        // return interpolation over momentum binning
        if(mom < fgTrackMomentum[0]) return lq1;
        else if(mom > fgTrackMomentum[kNMom-1]) return lq2;
        else return lq1 + (lq2 - lq1)*(mom - mom1)/(mom2 - mom1);

}

//_________________________________________________________________________
void AliTRDCalPIDLQ::Init()
{
  //
  // Initialization
  //

  fModel = new TObjArray(AliPID::kSPECIES  * kNMom);
  fModel -> SetOwner();

}

//_________________________________________________________________________
Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t /*ii*/) const
{  
  // returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)

  return spec * AliTRDCalPID::kNMom + mom;
}
Commit	Line	Data
44dbae42	1	/**************************************************************************
e7d6a389	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	**************************************************************************/
44dbae42	15
	16	/* $Id$ */
	17
	18	//////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Container for the distributions of dE/dx and the time bin of the //
	21	// max. cluster for electrons and pions //
	22	// //
	23	// Authors: //
	24	// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> //
	25	// Alex Bercuci <a.bercuci@gsi.de> //
	26	// //
	27	//////////////////////////////////////////////////////////////////////////
	28
44dbae42	29	#include <TH2F.h>
	30	#include <TFile.h>
	31	#include <TROOT.h>
	32
	33	#include "AliLog.h"
	34	#include "AliPID.h"
	35
e7d6a389	36	#include "../../STAT/TKDPDF.h"
44dbae42	37	#include "AliTRDCalPIDLQ.h"
	38	#include "AliTRDcalibDB.h"
	39
	40	ClassImp(AliTRDCalPIDLQ)
	41
2e32a5ae	42	Float_t AliTRDCalPIDLQ::fgTrackSegLength[kNLength] = { 3.7, 3.9, 4.2, 5.0 };
44dbae42	43
	44	//_________________________________________________________________________
	45	AliTRDCalPIDLQ::AliTRDCalPIDLQ()
	46	:AliTRDCalPID("pid", "LQ PID references for TRD")
	47	{
	48	//
	49	// The Default constructor
	50	//
	51
1fcc78b7	52	//Init();
44dbae42	53
	54	}
	55
	56	//_________________________________________________________________________
	57	AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t name, const Text_t title)
	58	:AliTRDCalPID(name,title)
	59	{
	60	//
	61	// The main constructor
	62	//
	63
	64	Init();
	65
	66	}
	67
	68	//_________________________________________________________________________
	69	AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
	70	{
	71	//
	72	// Destructor
	73	//
	74
	75	}
	76
e7d6a389	77	//_________________________________________________________________________
	78	Bool_t AliTRDCalPIDLQ::LoadPDF(TDirectoryFile *d)
	79	{
	80	// Read histograms
	81	TKDPDF *pdf = 0x0;
	82	for (Int_t is=0; is < AliPID::kSPECIES; is++){
	83	for (Int_t ip = 0; ip < kNMom; ip++){
	84	if(!(pdf = (TKDPDF*)d->Get(Form("%s[%d]", AliPID::ParticleShortName(is), ip)))){
	85	AliWarning(Form("Reference for %s[%d] missing.", AliPID::ParticleShortName(is), ip));
	86	continue;
	87	}
	88	fModel->AddAt(pdf->Clone(), GetModelID(ip, is, 0));
	89	}
	90	}
4f9b0c90	91	return kTRUE;
e7d6a389	92	}
e7d6a389	93	//
44dbae42	94	//_________________________________________________________________________
	95	Bool_t AliTRDCalPIDLQ::LoadReferences(Char_t *refFile)
	96	{
	97	//
	98	// Read the TRD dEdx histograms.
	99	//
	100
e7d6a389	101	Int_t nTimeBins = 22;
	102	// Get number of time bins from CDB
	103	AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
	104	if(!calibration){
	105	AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
	106	}else{
	107	if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
	108	else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
	109	}
	110
	111
44dbae42	112	// Read histogram Root file
	113	TFile *histFile = TFile::Open(refFile, "READ");
	114	if (!histFile \|\| !histFile->IsOpen()) {
	115	AliError(Form("Opening TRD histgram file %s failed", refFile));
	116	return kFALSE;
	117	}
	118	gROOT->cd();
	119
	120	// Read histograms
	121	for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
	122	for (Int_t imom = 0; imom < kNMom; imom++){
e7d6a389	123	TH2D* hist = (TH2D)histFile->Get(Form("h2dEdx%s%d", fPartSymb[iparticle], imom/, ilength*/))->Clone();
	124	hist->Scale(1./hist->Integral());
	125	fModel->AddAt(hist, GetModelID(imom, iparticle, 0));
44dbae42	126
	127	// if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
	128	//
	129	// TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fPartSymb[iparticle], imom))->Clone();
	130	// if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fPartSymb[iparticle], imom, nTimeBins));
	131	// ht->Scale(1./ht->Integral());
	132	// fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
e7d6a389	133	}
44dbae42	134	}
	135
	136	histFile->Close();
	137	delete histFile;
	138
	139	// Number of bins and bin size
	140	//TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
	141	//fNbins = hist->GetNbinsX();
	142	//fBinSize = hist->GetBinWidth(1);
	143
	144	return kTRUE;
	145
	146
	147	}
	148
	149	//_________________________________________________________________________
	150	TObject* AliTRDCalPIDLQ::GetModel(Int_t ip, Int_t iType, Int_t iplane) const // iType not needed
	151	{
	152	//
	153	// Returns one selected dEdx histogram
	154	//
	155
	156	if (iType < 0 \|\| iType >= AliPID::kSPECIES) return 0x0;
e7d6a389	157	if(ip<0 \|\| ip>= kNMom ) return 0x0;
44dbae42	158
2e32a5ae	159	AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fPartName[iType], fgTrackMomentum[ip]));
44dbae42	160
	161	return fModel->At(GetModelID(ip, iType, iplane));
	162	}
	163
	164	//_________________________________________________________________________
1fcc78b7	165	Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom
	166	, const Float_t * const dedx
	167	, Float_t length, Int_t iplane) const
44dbae42	168	{
44dbae42	169	//
e7d6a389	170	// Core function of AliTRDCalPID class for calculating the
	171	// likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
	172	// given momentum "mom" and a given dE/dx (slice "dedx") yield per
	173	// layer
44dbae42	174	//
44dbae42	175
e7d6a389	176	if (spec < 0 \|\| spec >= AliPID::kSPECIES) return 0.;
	177
	178	//Double_t dedx = dedx1/fMeanChargeRatio;
	179
	180	// find the interval in momentum and track segment length which applies for this data
	181	Int_t ilength = 1;
2e32a5ae	182	while(ilength<kNLength-1 && length>fgTrackSegLength[ilength]){
e7d6a389	183	ilength++;
	184	}
	185	Int_t imom = 1;
2e32a5ae	186	while(imom<kNMom-1 && mom>fgTrackMomentum[imom]) imom++;
44dbae42	187
e7d6a389	188	Int_t nbinsx, nbinsy;
	189	TAxis ax = 0x0, ay = 0x0;
	190	Double_t lq1, lq2;
	191	Double_t mom1 = fgTrackMomentum[imom-1], mom2 = fgTrackMomentum[imom];
	192	TH2 *hist = 0x0;
	193	if(!(hist = (TH2D*)fModel->At(GetModelID(imom-1, spec, iplane)))){
	194	AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
	195	AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom-1, spec, iplane)));
	196	return 0.;
	197	}
	198	ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
	199	ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
	200	Float_t x = dedx[0]+dedx[1], y = dedx[2];
44dbae42	201	Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
e7d6a389	202	Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
	203	if(kX)
	204	if(kY) lq1 = hist->GetBinContent( hist->FindBin(x, y));
44dbae42	205	//fEstimator->Estimate2D2(hist, x, y);
e7d6a389	206	else lq1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
	207	else
	208	if(kY) lq1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
	209	else lq1 = hist->GetBinContent(nbinsx, nbinsy);
	210
	211
	212	if(!(hist = (TH2D*)fModel->At(GetModelID(imom, spec, iplane)))){
	213	AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
	214	AliError(Form("EHistogram id %d not found in DB.", GetModelID(imom, spec, iplane)));
	215	return lq1;
	216	}
	217	if(kX)
	218	if(kY) lq2 = hist->GetBinContent( hist->FindBin(x, y));
44dbae42	219	//fEstimator->Estimate2D2(hist, x, y);
e7d6a389	220	else lq2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
	221	else
	222	if(kY) lq2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
	223	else lq2 = hist->GetBinContent(nbinsx, nbinsy);
	224
44dbae42	225	// return interpolation over momentum binning
2e32a5ae	226	if(mom < fgTrackMomentum[0]) return lq1;
	227	else if(mom > fgTrackMomentum[kNMom-1]) return lq2;
	228	else return lq1 + (lq2 - lq1)*(mom - mom1)/(mom2 - mom1);
44dbae42	229
	230	}
	231
	232	//_________________________________________________________________________
	233	void AliTRDCalPIDLQ::Init()
	234	{
	235	//
	236	// Initialization
	237	//
	238
	239	fModel = new TObjArray(AliPID::kSPECIES * kNMom);
	240	fModel -> SetOwner();
	241
	242	}
	243
	244	//_________________________________________________________________________
2e32a5ae	245	Int_t AliTRDCalPIDLQ::GetModelID(Int_t mom, Int_t spec, Int_t /ii/) const
44dbae42	246	{
	247	// returns the ID of the LQ distribution (55 Histos, ID from 1 to 55)
	248
e7d6a389	249	return spec * AliTRDCalPID::kNMom + mom;
44dbae42	250	}