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

/*
$Log$
Revision 1.2.8.2  2002/07/24 10:09:31  alibrary
Updating VirtualMC

Revision 1.2.8.1  2002/06/10 15:28:58  hristov
Merged with v3-08-02

Revision 1.4  2002/03/28 15:44:55  cblume
Remove const from GetIndex()

Revision 1.3  2002/03/28 14:59:07  cblume
Coding conventions

Revision 1.4  2002/03/28 15:44:55  cblume
Remove const from GetIndex()

Revision 1.3  2002/03/28 14:59:07  cblume
Coding conventions

Revision 1.2  2001/11/07 11:04:22  hristov
Minor corrections needed on Sun (arrays with undefined size created by new, inline decration removed when the body was hot in the header file)

Revision 1.1  2001/11/06 17:19:41  cblume
Add detailed geometry and simple simulator

*/

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//   The TRD particle identification class                                   //
//                                                                           //
//   Its main purposes are:                                                  //
//      - Creation and bookkeeping of the propability distributions          //
//      - Assignment of a e/pi propability to a given track based on         //
//        the LQ method                                                      //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <stdlib.h>
#include <math.h>

#include <TROOT.h>
#include <TH1.h>
#include <TObjArray.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>

#include "AliRun.h"
#include "AliTRD.h"
#include "AliTRDpidLQ.h"
#include "AliTRDcluster.h"
#include "AliTRDtrack.h"
#include "AliTRDtracker.h"
#include "AliTRDgeometry.h"

ClassImp(AliTRDpidLQ)

//_____________________________________________________________________________
AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
{
  //
  // AliTRDpidLQ default constructor
  // 

  fNMom   = 0;
  fMinMom = 0;
  fMaxMom = 0;
  fWidMom = 0;

  fNLh    = 0;
  fMinLh  = 0;
  fMaxLh  = 0;

  fHist   = NULL;

  fChargeMin   = 0;
  fNClusterMin = 0;

}

//_____________________________________________________________________________
AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
            :AliTRDpid(name,title)
{
  //
  // AliTRDpidLQ constructor
  //

  fNMom   = 0;
  fMinMom = 0;
  fMaxMom = 0;
  fWidMom = 0;

  fNLh    = 0;
  fMinLh  = 0;
  fMaxLh  = 0;

  fHist   = NULL;

  Init();

}

//_____________________________________________________________________________
AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p)
{
  //
  // AliTRDpidLQ copy constructor
  //

  ((AliTRDpidLQ &) p).Copy(*this);

}

//_____________________________________________________________________________
AliTRDpidLQ::~AliTRDpidLQ()
{
  //
  // AliTRDpidLQ destructor
  //

  if (fHist) {
    fHist->Delete();
    delete fHist;
  }

}

//_____________________________________________________________________________
AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
{
  //
  // Assignment operator
  //

  if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
  return *this;

}

//_____________________________________________________________________________
void AliTRDpidLQ::Copy(TObject &p)
{
  //
  // Copy function
  //

  fHist->Copy(*((AliTRDpidLQ &) p).fHist);

  ((AliTRDpidLQ &) p).fNMom        = fNMom;
  ((AliTRDpidLQ &) p).fMinMom      = fMinMom;
  ((AliTRDpidLQ &) p).fMaxMom      = fMaxMom;
  ((AliTRDpidLQ &) p).fWidMom      = fWidMom;
  ((AliTRDpidLQ &) p).fChargeMin   = fChargeMin;
  ((AliTRDpidLQ &) p).fNClusterMin = fNClusterMin;
  ((AliTRDpidLQ &) p).fNLh         = fNLh;
  ((AliTRDpidLQ &) p).fMinLh       = fMinLh;
  ((AliTRDpidLQ &) p).fMaxLh       = fMaxLh;

  AliTRDpid::Copy(p);

}

//_____________________________________________________________________________
Bool_t AliTRDpidLQ::Init()
{
  //
  // Initializes the PID object 
  //

  fChargeMin   = 0.0;
  fNClusterMin = 10;

  fNLh         = 100;
  fMinLh       =   0.0;
  fMaxLh       = 500.0;

  return kTRUE;

}

//_____________________________________________________________________________
Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
{
  //
  // Assigns the e / pi likelihood to a given track
  //

  const Int_t kNpla = AliTRDgeometry::Nplan();
  Float_t * charge = new Float_t[kNpla];
  Int_t   * nCluster = new Int_t[kNpla];

  Float_t  lhPi = 0;
  Float_t  lhEl = 0;

  Double_t pPi  = 1;
  Double_t pEl  = 1;
  Double_t pSum = 0;

  Int_t   indexEl;
  Int_t   indexPi;
  TH1F   *hTmpEl;
  TH1F   *hTmpPi;

  t->SetLikelihoodElectron(-1.);
  if (isnan(t->GetP())) return kFALSE;
  Float_t mom = t->GetP();

  // Calculate the total charge in each plane
  if (!SumCharge(t,charge,nCluster)) return kFALSE;

  indexEl = GetIndex(mom,kElectron);
  indexPi = GetIndex(mom,kPion);
  if ((indexEl > -1) && (indexPi > -1)) {
    hTmpEl = (TH1F *) fHist->UncheckedAt(indexEl);
    hTmpPi = (TH1F *) fHist->UncheckedAt(indexPi);
    for (Int_t ipla = 0; ipla < kNpla; ipla++) {
      Float_t chargePlane   = charge[ipla];
      Int_t   nClusterPlane = nCluster[ipla];
      if ((chargePlane   >   fChargeMin) &&
          (nClusterPlane > fNClusterMin)){
        Float_t chargeNorm = chargePlane / ((Float_t) nClusterPlane);
        if (chargeNorm < fMaxLh) {
          Int_t   ibinEl     = hTmpEl->FindBin(chargeNorm);
          Float_t pElPlane   = hTmpEl->GetBinContent(ibinEl);
          if (pElPlane > 0) { 
            pEl  = pEl * pElPlane;
          } 
          Int_t   ibinPi     = hTmpPi->FindBin(chargeNorm);
          Float_t pPiPlane   = hTmpPi->GetBinContent(ibinPi);
          if (pPiPlane > 0) { 
            pPi  = pPi * pPiPlane;
          } 
//           printf("          ipla = %d, nCluster = %d, charge = %f\n"
//                  ,ipla,nClusterPlane,chargeNorm);
//           printf("electron: pElPlane = %f, ibinEl = %d, pEl = %f\n"
// 	         ,pElPlane,ibinEl,pEl);
//           printf("    pion: pPiPlane = %f, ibinPi = %d, pPi = %f\n"
// 	         ,pPiPlane,ibinPi,pPi);
	}
      }
    }
  }
  else {
    delete [] charge;
    delete [] nCluster;
    return kTRUE;
  }

  pSum = pEl + pPi;
  if (pSum <= 0) {
    delete [] charge;
    delete [] nCluster;
    return kFALSE;
  }
  lhEl = pEl / pSum;
  lhPi = pPi / pSum;

//   printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);

  // Assign the likelihoods 
  t->SetLikelihoodElectron(lhEl);

  delete [] charge;
  delete [] nCluster;
  return kTRUE;  

}

//_____________________________________________________________________________
Bool_t AliTRDpidLQ::CreateHistograms(const Int_t   nmom
                                   , const Float_t minmom
                                   , const Float_t maxmom)
{
  //
  // Creates the histograms
  //

  Int_t imom;
  Int_t ipid;

  fNMom   = nmom;
  fMinMom = minmom;
  fMaxMom = maxmom;
  fWidMom = (maxmom - minmom) / ((Float_t) nmom);

  fHist = new TObjArray(kNpid * nmom);
  for (imom = 0; imom < nmom;  imom++) {
    for (ipid = 0; ipid < kNpid; ipid++) {

      Int_t index = GetIndex(imom,ipid);
      Char_t name[10];
      Char_t title[80];
      sprintf(name ,"hQ%03d",index);
      if (ipid == kElectron) {
        sprintf(title,"Q electron p-bin %03d",imom);
      }
      else {
        sprintf(title,"Q pion p-bin %03d",imom);
      }
      TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
      fHist->AddAt(hTmp,index);

    }
  }

  return kTRUE;

}

// //_____________________________________________________________________________
// Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
// {
//   //
//   // Fills the energy loss distributions with track <t>.
//   //

//   Bool_t status = kTRUE;

//   if (isnan(t->GetP())) return kFALSE;

//   Float_t        mom     = t->GetP();
//   Int_t          ipid    = MCpid(t);
//   TH1F          *hTmp    = NULL;
//   AliTRDcluster *cluster = NULL;

//   Int_t index = GetIndex(mom,ipid);
//   if (index > -1) {
//     hTmp = (TH1F *) fHist->UncheckedAt(index);
//     // Loop through all clusters associated to this track
//     Int_t nClus = t->GetNclusters();
//     for (Int_t iClus = 0; iClus < nClus; iClus++) {
//       // Get a cluster
//       Int_t idxClus = t->GetClusterIndex(iClus);
//       if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
//         status = kFALSE;
//         break;
//       } 
//       hTmp->Fill(cluster->GetQ());
//     }
//   }  

//   return status;

// }

//_____________________________________________________________________________
Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
{
  //
  // Fills the energy loss distributions with track <t>.
  //

  const Int_t kNpla = AliTRDgeometry::Nplan();

  if (isnan(t->GetP())) return kFALSE;

  Float_t * charge = new Float_t[kNpla];
  Int_t   * nCluster = new Int_t[kNpla];
  Float_t mom  = t->GetP();
  Int_t   ipid = MCpid(t);
  TH1F   *hTmp = NULL;

  if (!SumCharge(t,charge,nCluster)) {
    delete [] charge;
    delete [] nCluster;
    return kFALSE;
  }

  Int_t index = GetIndex(mom,ipid);
  if (index > -1) {
    hTmp = (TH1F *) fHist->UncheckedAt(index);
    for (Int_t ipla = 0; ipla < kNpla; ipla++) {
      if ((charge[ipla]   >   fChargeMin) &&
          (nCluster[ipla] > fNClusterMin)){
        hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
      }
    }
  }  

  delete [] charge;
  delete [] nCluster;
  return kTRUE;

}

//_____________________________________________________________________________
Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
{
  //
  // Returns the histogram index
  //

  if (isnan(t->GetP())) return -1;
  Float_t mom  = t->GetP();
  Int_t   ipid = MCpid(t);

  return GetIndex(mom,ipid);

}

//_____________________________________________________________________________
Int_t AliTRDpidLQ::GetIndex(const Float_t mom, const Int_t ipid)
{
  //
  // Returns the histogram index
  //

  if ((mom < fMinMom) || (mom > fMaxMom))  return -1;
  Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
  return GetIndex(imom,ipid);

}

//_____________________________________________________________________________
Int_t AliTRDpidLQ::GetIndex(const Int_t imom, const Int_t ipid)
{
  //
  // Returns the histogram index
  //

  if ((ipid < 0) || (ipid >= kNpid)) return -1;
  return imom * kNpid + ipid; 

}
Commit	Line	Data
16bf9884	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	/*
	17	$Log$
b9d0a01d	18	Revision 1.2.8.2 2002/07/24 10:09:31 alibrary
	19	Updating VirtualMC
	20
	21	Revision 1.2.8.1 2002/06/10 15:28:58 hristov
	22	Merged with v3-08-02
	23
	24	Revision 1.4 2002/03/28 15:44:55 cblume
	25	Remove const from GetIndex()
	26
	27	Revision 1.3 2002/03/28 14:59:07 cblume
	28	Coding conventions
	29
	30	Revision 1.4 2002/03/28 15:44:55 cblume
	31	Remove const from GetIndex()
	32
95fc95f8	33	Revision 1.3 2002/03/28 14:59:07 cblume
	34	Coding conventions
	35
0a29d0f1	36	Revision 1.2 2001/11/07 11:04:22 hristov
	37	Minor corrections needed on Sun (arrays with undefined size created by new, inline decration removed when the body was hot in the header file)
	38
a3d851c2	39	Revision 1.1 2001/11/06 17:19:41 cblume
	40	Add detailed geometry and simple simulator
	41
16bf9884	42	*/
	43
	44	///////////////////////////////////////////////////////////////////////////////
	45	// //
	46	// The TRD particle identification class //
	47	// //
	48	// Its main purposes are: //
	49	// - Creation and bookkeeping of the propability distributions //
	50	// - Assignment of a e/pi propability to a given track based on //
	51	// the LQ method //
	52	// //
	53	///////////////////////////////////////////////////////////////////////////////
	54
	55	#include <stdlib.h>
	56	#include <math.h>
	57
	58	#include <TROOT.h>
	59	#include <TH1.h>
	60	#include <TObjArray.h>
	61	#include <TTree.h>
	62	#include <TFile.h>
	63	#include <TParticle.h>
	64
	65	#include "AliRun.h"
	66	#include "AliTRD.h"
	67	#include "AliTRDpidLQ.h"
	68	#include "AliTRDcluster.h"
	69	#include "AliTRDtrack.h"
	70	#include "AliTRDtracker.h"
	71	#include "AliTRDgeometry.h"
	72
	73	ClassImp(AliTRDpidLQ)
	74
	75	//_____________________________________________________________________________
	76	AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
	77	{
	78	//
	79	// AliTRDpidLQ default constructor
	80	//
	81
	82	fNMom = 0;
	83	fMinMom = 0;
	84	fMaxMom = 0;
	85	fWidMom = 0;
	86
	87	fNLh = 0;
	88	fMinLh = 0;
	89	fMaxLh = 0;
	90
	91	fHist = NULL;
	92
	93	fChargeMin = 0;
	94	fNClusterMin = 0;
	95
	96	}
	97
	98	//_____________________________________________________________________________
	99	AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
	100	:AliTRDpid(name,title)
	101	{
	102	//
	103	// AliTRDpidLQ constructor
	104	//
	105
106	fNMom = 0;
107	fMinMom = 0;
108	fMaxMom = 0;
109	fWidMom = 0;
110
111	fNLh = 0;
112	fMinLh = 0;
113	fMaxLh = 0;
114
115	fHist = NULL;
116
117	Init();
118
119	}
120
121	//_____________________________________________________________________________
122	AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p)
123	{
124	//
125	// AliTRDpidLQ copy constructor
126	//
127
128	((AliTRDpidLQ &) p).Copy(*this);
129
130	}
131
132	//_____________________________________________________________________________
133	AliTRDpidLQ::~AliTRDpidLQ()
134	{
135	//
136	// AliTRDpidLQ destructor
137	//
138
139	if (fHist) {
140	fHist->Delete();
141	delete fHist;
142	}
143
144	}
145
146	//_____________________________________________________________________________
147	AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
148	{
149	//
150	// Assignment operator
151	//
152
153	if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
154	return *this;
155
156	}
157
158	//_____________________________________________________________________________
159	void AliTRDpidLQ::Copy(TObject &p)
160	{
161	//
162	// Copy function
163	//
164
165	fHist->Copy(*((AliTRDpidLQ &) p).fHist);
166
167	((AliTRDpidLQ &) p).fNMom = fNMom;
168	((AliTRDpidLQ &) p).fMinMom = fMinMom;
169	((AliTRDpidLQ &) p).fMaxMom = fMaxMom;
170	((AliTRDpidLQ &) p).fWidMom = fWidMom;
171	((AliTRDpidLQ &) p).fChargeMin = fChargeMin;
172	((AliTRDpidLQ &) p).fNClusterMin = fNClusterMin;
173	((AliTRDpidLQ &) p).fNLh = fNLh;
174	((AliTRDpidLQ &) p).fMinLh = fMinLh;
175	((AliTRDpidLQ &) p).fMaxLh = fMaxLh;
176
177	AliTRDpid::Copy(p);
178
179	}
180
181	//_____________________________________________________________________________
182	Bool_t AliTRDpidLQ::Init()
183	{
184	//
185	// Initializes the PID object
186	//
187
188	fChargeMin = 0.0;
189	fNClusterMin = 10;
190
191	fNLh = 100;
192	fMinLh = 0.0;
193	fMaxLh = 500.0;
194
195	return kTRUE;
196
197	}
198
199	//_____________________________________________________________________________
200	Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
201	{
202	//
203	// Assigns the e / pi likelihood to a given track
204	//
205
206	const Int_t kNpla = AliTRDgeometry::Nplan();
a3d851c2	207	Float_t * charge = new Float_t[kNpla];
a3d851c2	208	Int_t * nCluster = new Int_t[kNpla];
16bf9884	209
	210	Float_t lhPi = 0;
	211	Float_t lhEl = 0;
	212
	213	Double_t pPi = 1;
	214	Double_t pEl = 1;
	215	Double_t pSum = 0;
	216
	217	Int_t indexEl;
	218	Int_t indexPi;
	219	TH1F *hTmpEl;
	220	TH1F *hTmpPi;
	221
	222	t->SetLikelihoodElectron(-1.);
	223	if (isnan(t->GetP())) return kFALSE;
	224	Float_t mom = t->GetP();
	225
	226	// Calculate the total charge in each plane
	227	if (!SumCharge(t,charge,nCluster)) return kFALSE;
	228
	229	indexEl = GetIndex(mom,kElectron);
	230	indexPi = GetIndex(mom,kPion);
	231	if ((indexEl > -1) && (indexPi > -1)) {
	232	hTmpEl = (TH1F *) fHist->UncheckedAt(indexEl);
	233	hTmpPi = (TH1F *) fHist->UncheckedAt(indexPi);
	234	for (Int_t ipla = 0; ipla < kNpla; ipla++) {
	235	Float_t chargePlane = charge[ipla];
	236	Int_t nClusterPlane = nCluster[ipla];
	237	if ((chargePlane > fChargeMin) &&
	238	(nClusterPlane > fNClusterMin)){
	239	Float_t chargeNorm = chargePlane / ((Float_t) nClusterPlane);
	240	if (chargeNorm < fMaxLh) {
	241	Int_t ibinEl = hTmpEl->FindBin(chargeNorm);
	242	Float_t pElPlane = hTmpEl->GetBinContent(ibinEl);
	243	if (pElPlane > 0) {
	244	pEl = pEl * pElPlane;
	245	}
	246	Int_t ibinPi = hTmpPi->FindBin(chargeNorm);
	247	Float_t pPiPlane = hTmpPi->GetBinContent(ibinPi);
	248	if (pPiPlane > 0) {
	249	pPi = pPi * pPiPlane;
	250	}
	251	// printf(" ipla = %d, nCluster = %d, charge = %f\n"
	252	// ,ipla,nClusterPlane,chargeNorm);
	253	// printf("electron: pElPlane = %f, ibinEl = %d, pEl = %f\n"
	254	// ,pElPlane,ibinEl,pEl);
	255	// printf(" pion: pPiPlane = %f, ibinPi = %d, pPi = %f\n"
	256	// ,pPiPlane,ibinPi,pPi);
	257	}
	258	}
	259	}
	260	}
	261	else {
a3d851c2	262	delete [] charge;
a3d851c2	263	delete [] nCluster;
16bf9884	264	return kTRUE;
	265	}
	266
	267	pSum = pEl + pPi;
a3d851c2	268	if (pSum <= 0) {
	269	delete [] charge;
	270	delete [] nCluster;
	271	return kFALSE;
	272	}
16bf9884	273	lhEl = pEl / pSum;
	274	lhPi = pPi / pSum;
	275
	276	// printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);
	277
	278	// Assign the likelihoods
	279	t->SetLikelihoodElectron(lhEl);
	280
a3d851c2	281	delete [] charge;
a3d851c2	282	delete [] nCluster;
16bf9884	283	return kTRUE;
	284
	285	}
	286
	287	//_____________________________________________________________________________
	288	Bool_t AliTRDpidLQ::CreateHistograms(const Int_t nmom
	289	, const Float_t minmom
	290	, const Float_t maxmom)
	291	{
	292	//
	293	// Creates the histograms
	294	//
	295
	296	Int_t imom;
	297	Int_t ipid;
	298
	299	fNMom = nmom;
	300	fMinMom = minmom;
	301	fMaxMom = maxmom;
	302	fWidMom = (maxmom - minmom) / ((Float_t) nmom);
	303
	304	fHist = new TObjArray(kNpid * nmom);
	305	for (imom = 0; imom < nmom; imom++) {
	306	for (ipid = 0; ipid < kNpid; ipid++) {
	307
	308	Int_t index = GetIndex(imom,ipid);
	309	Char_t name[10];
	310	Char_t title[80];
	311	sprintf(name ,"hQ%03d",index);
	312	if (ipid == kElectron) {
	313	sprintf(title,"Q electron p-bin %03d",imom);
	314	}
	315	else {
	316	sprintf(title,"Q pion p-bin %03d",imom);
	317	}
	318	TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
	319	fHist->AddAt(hTmp,index);
	320
	321	}
	322	}
	323
	324	return kTRUE;
	325
	326	}
	327
	328	// //_____________________________________________________________________________
	329	// Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
	330	// {
	331	// //
	332	// // Fills the energy loss distributions with track <t>.
	333	// //
	334
	335	// Bool_t status = kTRUE;
	336
	337	// if (isnan(t->GetP())) return kFALSE;
	338
	339	// Float_t mom = t->GetP();
	340	// Int_t ipid = MCpid(t);
	341	// TH1F *hTmp = NULL;
	342	// AliTRDcluster *cluster = NULL;
	343
	344	// Int_t index = GetIndex(mom,ipid);
	345	// if (index > -1) {
	346	// hTmp = (TH1F *) fHist->UncheckedAt(index);
347	// // Loop through all clusters associated to this track
348	// Int_t nClus = t->GetNclusters();
349	// for (Int_t iClus = 0; iClus < nClus; iClus++) {
350	// // Get a cluster
351	// Int_t idxClus = t->GetClusterIndex(iClus);
352	// if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
353	// status = kFALSE;
354	// break;
355	// }
356	// hTmp->Fill(cluster->GetQ());
357	// }
358	// }
359
360	// return status;
361
362	// }
363
364	//_____________________________________________________________________________
365	Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
366	{
367	//
368	// Fills the energy loss distributions with track <t>.
369	//
370
371	const Int_t kNpla = AliTRDgeometry::Nplan();
372
373	if (isnan(t->GetP())) return kFALSE;
374
a3d851c2	375	Float_t * charge = new Float_t[kNpla];
a3d851c2	376	Int_t * nCluster = new Int_t[kNpla];
16bf9884	377	Float_t mom = t->GetP();
	378	Int_t ipid = MCpid(t);
	379	TH1F *hTmp = NULL;
	380
a3d851c2	381	if (!SumCharge(t,charge,nCluster)) {
	382	delete [] charge;
	383	delete [] nCluster;
	384	return kFALSE;
	385	}
16bf9884	386
	387	Int_t index = GetIndex(mom,ipid);
	388	if (index > -1) {
	389	hTmp = (TH1F *) fHist->UncheckedAt(index);
	390	for (Int_t ipla = 0; ipla < kNpla; ipla++) {
	391	if ((charge[ipla] > fChargeMin) &&
	392	(nCluster[ipla] > fNClusterMin)){
	393	hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
	394	}
	395	}
	396	}
	397
a3d851c2	398	delete [] charge;
a3d851c2	399	delete [] nCluster;
16bf9884	400	return kTRUE;
	401
	402	}
	403
	404	//_____________________________________________________________________________
	405	Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
	406	{
	407	//
	408	// Returns the histogram index
	409	//
	410
	411	if (isnan(t->GetP())) return -1;
	412	Float_t mom = t->GetP();
	413	Int_t ipid = MCpid(t);
	414
	415	return GetIndex(mom,ipid);
	416
	417	}
	418
	419	//_____________________________________________________________________________
	420	Int_t AliTRDpidLQ::GetIndex(const Float_t mom, const Int_t ipid)
	421	{
	422	//
	423	// Returns the histogram index
	424	//
	425
	426	if ((mom < fMinMom) \|\| (mom > fMaxMom)) return -1;
	427	Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
	428	return GetIndex(imom,ipid);
	429
	430	}
	431
	432	//_____________________________________________________________________________
95fc95f8	433	Int_t AliTRDpidLQ::GetIndex(const Int_t imom, const Int_t ipid)
16bf9884	434	{
	435	//
	436	// Returns the histogram index
	437	//
	438
	439	if ((ipid < 0) \|\| (ipid >= kNpid)) return -1;
	440	return imom * kNpid + ipid;
	441
	442	}