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