1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.3 2002/03/28 14:59:07 cblume
21 Revision 1.2 2001/11/07 11:04:22 hristov
22 Minor corrections needed on Sun (arrays with undefined size created by new, inline decration removed when the body was hot in the header file)
24 Revision 1.1 2001/11/06 17:19:41 cblume
25 Add detailed geometry and simple simulator
29 ///////////////////////////////////////////////////////////////////////////////
31 // The TRD particle identification class //
33 // Its main purposes are: //
34 // - Creation and bookkeeping of the propability distributions //
35 // - Assignment of a e/pi propability to a given track based on //
38 ///////////////////////////////////////////////////////////////////////////////
45 #include <TObjArray.h>
48 #include <TParticle.h>
52 #include "AliTRDpidLQ.h"
53 #include "AliTRDcluster.h"
54 #include "AliTRDtrack.h"
55 #include "AliTRDtracker.h"
56 #include "AliTRDgeometry.h"
60 //_____________________________________________________________________________
61 AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
64 // AliTRDpidLQ default constructor
83 //_____________________________________________________________________________
84 AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
85 :AliTRDpid(name,title)
88 // AliTRDpidLQ constructor
106 //_____________________________________________________________________________
107 AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p)
110 // AliTRDpidLQ copy constructor
113 ((AliTRDpidLQ &) p).Copy(*this);
117 //_____________________________________________________________________________
118 AliTRDpidLQ::~AliTRDpidLQ()
121 // AliTRDpidLQ destructor
131 //_____________________________________________________________________________
132 AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
135 // Assignment operator
138 if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
143 //_____________________________________________________________________________
144 void AliTRDpidLQ::Copy(TObject &p)
150 fHist->Copy(*((AliTRDpidLQ &) p).fHist);
152 ((AliTRDpidLQ &) p).fNMom = fNMom;
153 ((AliTRDpidLQ &) p).fMinMom = fMinMom;
154 ((AliTRDpidLQ &) p).fMaxMom = fMaxMom;
155 ((AliTRDpidLQ &) p).fWidMom = fWidMom;
156 ((AliTRDpidLQ &) p).fChargeMin = fChargeMin;
157 ((AliTRDpidLQ &) p).fNClusterMin = fNClusterMin;
158 ((AliTRDpidLQ &) p).fNLh = fNLh;
159 ((AliTRDpidLQ &) p).fMinLh = fMinLh;
160 ((AliTRDpidLQ &) p).fMaxLh = fMaxLh;
166 //_____________________________________________________________________________
167 Bool_t AliTRDpidLQ::Init()
170 // Initializes the PID object
184 //_____________________________________________________________________________
185 Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
188 // Assigns the e / pi likelihood to a given track
191 const Int_t kNpla = AliTRDgeometry::Nplan();
192 Float_t * charge = new Float_t[kNpla];
193 Int_t * nCluster = new Int_t[kNpla];
207 t->SetLikelihoodElectron(-1.);
208 if (isnan(t->GetP())) return kFALSE;
209 Float_t mom = t->GetP();
211 // Calculate the total charge in each plane
212 if (!SumCharge(t,charge,nCluster)) return kFALSE;
214 indexEl = GetIndex(mom,kElectron);
215 indexPi = GetIndex(mom,kPion);
216 if ((indexEl > -1) && (indexPi > -1)) {
217 hTmpEl = (TH1F *) fHist->UncheckedAt(indexEl);
218 hTmpPi = (TH1F *) fHist->UncheckedAt(indexPi);
219 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
220 Float_t chargePlane = charge[ipla];
221 Int_t nClusterPlane = nCluster[ipla];
222 if ((chargePlane > fChargeMin) &&
223 (nClusterPlane > fNClusterMin)){
224 Float_t chargeNorm = chargePlane / ((Float_t) nClusterPlane);
225 if (chargeNorm < fMaxLh) {
226 Int_t ibinEl = hTmpEl->FindBin(chargeNorm);
227 Float_t pElPlane = hTmpEl->GetBinContent(ibinEl);
229 pEl = pEl * pElPlane;
231 Int_t ibinPi = hTmpPi->FindBin(chargeNorm);
232 Float_t pPiPlane = hTmpPi->GetBinContent(ibinPi);
234 pPi = pPi * pPiPlane;
236 // printf(" ipla = %d, nCluster = %d, charge = %f\n"
237 // ,ipla,nClusterPlane,chargeNorm);
238 // printf("electron: pElPlane = %f, ibinEl = %d, pEl = %f\n"
239 // ,pElPlane,ibinEl,pEl);
240 // printf(" pion: pPiPlane = %f, ibinPi = %d, pPi = %f\n"
241 // ,pPiPlane,ibinPi,pPi);
261 // printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);
263 // Assign the likelihoods
264 t->SetLikelihoodElectron(lhEl);
272 //_____________________________________________________________________________
273 Bool_t AliTRDpidLQ::CreateHistograms(const Int_t nmom
274 , const Float_t minmom
275 , const Float_t maxmom)
278 // Creates the histograms
287 fWidMom = (maxmom - minmom) / ((Float_t) nmom);
289 fHist = new TObjArray(kNpid * nmom);
290 for (imom = 0; imom < nmom; imom++) {
291 for (ipid = 0; ipid < kNpid; ipid++) {
293 Int_t index = GetIndex(imom,ipid);
296 sprintf(name ,"hQ%03d",index);
297 if (ipid == kElectron) {
298 sprintf(title,"Q electron p-bin %03d",imom);
301 sprintf(title,"Q pion p-bin %03d",imom);
303 TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
304 fHist->AddAt(hTmp,index);
313 // //_____________________________________________________________________________
314 // Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
317 // // Fills the energy loss distributions with track <t>.
320 // Bool_t status = kTRUE;
322 // if (isnan(t->GetP())) return kFALSE;
324 // Float_t mom = t->GetP();
325 // Int_t ipid = MCpid(t);
326 // TH1F *hTmp = NULL;
327 // AliTRDcluster *cluster = NULL;
329 // Int_t index = GetIndex(mom,ipid);
331 // hTmp = (TH1F *) fHist->UncheckedAt(index);
332 // // Loop through all clusters associated to this track
333 // Int_t nClus = t->GetNclusters();
334 // for (Int_t iClus = 0; iClus < nClus; iClus++) {
336 // Int_t idxClus = t->GetClusterIndex(iClus);
337 // if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
341 // hTmp->Fill(cluster->GetQ());
349 //_____________________________________________________________________________
350 Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
353 // Fills the energy loss distributions with track <t>.
356 const Int_t kNpla = AliTRDgeometry::Nplan();
358 if (isnan(t->GetP())) return kFALSE;
360 Float_t * charge = new Float_t[kNpla];
361 Int_t * nCluster = new Int_t[kNpla];
362 Float_t mom = t->GetP();
363 Int_t ipid = MCpid(t);
366 if (!SumCharge(t,charge,nCluster)) {
372 Int_t index = GetIndex(mom,ipid);
374 hTmp = (TH1F *) fHist->UncheckedAt(index);
375 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
376 if ((charge[ipla] > fChargeMin) &&
377 (nCluster[ipla] > fNClusterMin)){
378 hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
389 //_____________________________________________________________________________
390 Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
393 // Returns the histogram index
396 if (isnan(t->GetP())) return -1;
397 Float_t mom = t->GetP();
398 Int_t ipid = MCpid(t);
400 return GetIndex(mom,ipid);
404 //_____________________________________________________________________________
405 Int_t AliTRDpidLQ::GetIndex(const Float_t mom, const Int_t ipid)
408 // Returns the histogram index
411 if ((mom < fMinMom) || (mom > fMaxMom)) return -1;
412 Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
413 return GetIndex(imom,ipid);
417 //_____________________________________________________________________________
418 Int_t AliTRDpidLQ::GetIndex(const Int_t imom, const Int_t ipid)
421 // Returns the histogram index
424 if ((ipid < 0) || (ipid >= kNpid)) return -1;
425 return imom * kNpid + ipid;