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.1 2001/11/06 17:19:41 cblume
19 Add detailed geometry and simple simulator
23 ///////////////////////////////////////////////////////////////////////////////
25 // The TRD particle identification class //
27 // Its main purposes are: //
28 // - Creation and bookkeeping of the propability distributions //
29 // - Assignment of a e/pi propability to a given track based on //
32 ///////////////////////////////////////////////////////////////////////////////
39 #include <TObjArray.h>
42 #include <TParticle.h>
46 #include "AliTRDpidLQ.h"
47 #include "AliTRDcluster.h"
48 #include "AliTRDtrack.h"
49 #include "AliTRDtracker.h"
50 #include "AliTRDgeometry.h"
54 //_____________________________________________________________________________
55 AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
58 // AliTRDpidLQ default constructor
77 //_____________________________________________________________________________
78 AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
79 :AliTRDpid(name,title)
82 // AliTRDpidLQ constructor
100 //_____________________________________________________________________________
101 AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p)
104 // AliTRDpidLQ copy constructor
107 ((AliTRDpidLQ &) p).Copy(*this);
111 //_____________________________________________________________________________
112 AliTRDpidLQ::~AliTRDpidLQ()
115 // AliTRDpidLQ destructor
125 //_____________________________________________________________________________
126 AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
129 // Assignment operator
132 if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
137 //_____________________________________________________________________________
138 void AliTRDpidLQ::Copy(TObject &p)
144 fHist->Copy(*((AliTRDpidLQ &) p).fHist);
146 ((AliTRDpidLQ &) p).fNMom = fNMom;
147 ((AliTRDpidLQ &) p).fMinMom = fMinMom;
148 ((AliTRDpidLQ &) p).fMaxMom = fMaxMom;
149 ((AliTRDpidLQ &) p).fWidMom = fWidMom;
150 ((AliTRDpidLQ &) p).fChargeMin = fChargeMin;
151 ((AliTRDpidLQ &) p).fNClusterMin = fNClusterMin;
152 ((AliTRDpidLQ &) p).fNLh = fNLh;
153 ((AliTRDpidLQ &) p).fMinLh = fMinLh;
154 ((AliTRDpidLQ &) p).fMaxLh = fMaxLh;
160 //_____________________________________________________________________________
161 Bool_t AliTRDpidLQ::Init()
164 // Initializes the PID object
178 //_____________________________________________________________________________
179 Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
182 // Assigns the e / pi likelihood to a given track
185 const Int_t kNpla = AliTRDgeometry::Nplan();
186 Float_t * charge = new Float_t[kNpla];
187 Int_t * nCluster = new Int_t[kNpla];
201 t->SetLikelihoodElectron(-1.);
202 if (isnan(t->GetP())) return kFALSE;
203 Float_t mom = t->GetP();
205 // Calculate the total charge in each plane
206 if (!SumCharge(t,charge,nCluster)) return kFALSE;
208 indexEl = GetIndex(mom,kElectron);
209 indexPi = GetIndex(mom,kPion);
210 if ((indexEl > -1) && (indexPi > -1)) {
211 hTmpEl = (TH1F *) fHist->UncheckedAt(indexEl);
212 hTmpPi = (TH1F *) fHist->UncheckedAt(indexPi);
213 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
214 Float_t chargePlane = charge[ipla];
215 Int_t nClusterPlane = nCluster[ipla];
216 if ((chargePlane > fChargeMin) &&
217 (nClusterPlane > fNClusterMin)){
218 Float_t chargeNorm = chargePlane / ((Float_t) nClusterPlane);
219 if (chargeNorm < fMaxLh) {
220 Int_t ibinEl = hTmpEl->FindBin(chargeNorm);
221 Float_t pElPlane = hTmpEl->GetBinContent(ibinEl);
223 pEl = pEl * pElPlane;
225 Int_t ibinPi = hTmpPi->FindBin(chargeNorm);
226 Float_t pPiPlane = hTmpPi->GetBinContent(ibinPi);
228 pPi = pPi * pPiPlane;
230 // printf(" ipla = %d, nCluster = %d, charge = %f\n"
231 // ,ipla,nClusterPlane,chargeNorm);
232 // printf("electron: pElPlane = %f, ibinEl = %d, pEl = %f\n"
233 // ,pElPlane,ibinEl,pEl);
234 // printf(" pion: pPiPlane = %f, ibinPi = %d, pPi = %f\n"
235 // ,pPiPlane,ibinPi,pPi);
255 // printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);
257 // Assign the likelihoods
258 t->SetLikelihoodElectron(lhEl);
266 //_____________________________________________________________________________
267 Bool_t AliTRDpidLQ::CreateHistograms(const Int_t nmom
268 , const Float_t minmom
269 , const Float_t maxmom)
272 // Creates the histograms
281 fWidMom = (maxmom - minmom) / ((Float_t) nmom);
283 fHist = new TObjArray(kNpid * nmom);
284 for (imom = 0; imom < nmom; imom++) {
285 for (ipid = 0; ipid < kNpid; ipid++) {
287 Int_t index = GetIndex(imom,ipid);
290 sprintf(name ,"hQ%03d",index);
291 if (ipid == kElectron) {
292 sprintf(title,"Q electron p-bin %03d",imom);
295 sprintf(title,"Q pion p-bin %03d",imom);
297 TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
298 fHist->AddAt(hTmp,index);
307 // //_____________________________________________________________________________
308 // Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
311 // // Fills the energy loss distributions with track <t>.
314 // Bool_t status = kTRUE;
316 // if (isnan(t->GetP())) return kFALSE;
318 // Float_t mom = t->GetP();
319 // Int_t ipid = MCpid(t);
320 // TH1F *hTmp = NULL;
321 // AliTRDcluster *cluster = NULL;
323 // Int_t index = GetIndex(mom,ipid);
325 // hTmp = (TH1F *) fHist->UncheckedAt(index);
326 // // Loop through all clusters associated to this track
327 // Int_t nClus = t->GetNclusters();
328 // for (Int_t iClus = 0; iClus < nClus; iClus++) {
330 // Int_t idxClus = t->GetClusterIndex(iClus);
331 // if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
335 // hTmp->Fill(cluster->GetQ());
343 //_____________________________________________________________________________
344 Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
347 // Fills the energy loss distributions with track <t>.
350 const Int_t kNpla = AliTRDgeometry::Nplan();
352 if (isnan(t->GetP())) return kFALSE;
354 Float_t * charge = new Float_t[kNpla];
355 Int_t * nCluster = new Int_t[kNpla];
356 Float_t mom = t->GetP();
357 Int_t ipid = MCpid(t);
360 if (!SumCharge(t,charge,nCluster)) {
366 Int_t index = GetIndex(mom,ipid);
368 hTmp = (TH1F *) fHist->UncheckedAt(index);
369 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
370 if ((charge[ipla] > fChargeMin) &&
371 (nCluster[ipla] > fNClusterMin)){
372 hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
383 //_____________________________________________________________________________
384 Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
387 // Returns the histogram index
390 if (isnan(t->GetP())) return -1;
391 Float_t mom = t->GetP();
392 Int_t ipid = MCpid(t);
394 return GetIndex(mom,ipid);
398 //_____________________________________________________________________________
399 Int_t AliTRDpidLQ::GetIndex(const Float_t mom, const Int_t ipid)
402 // Returns the histogram index
405 if ((mom < fMinMom) || (mom > fMaxMom)) return -1;
406 Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
407 return GetIndex(imom,ipid);
411 //_____________________________________________________________________________
412 Int_t AliTRDpidLQ::GetIndex(const Int_t imom, const Int_t ipid)
415 // Returns the histogram index
418 if ((ipid < 0) || (ipid >= kNpid)) return -1;
419 return imom * kNpid + ipid;