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

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// AliHMPIDPid                                                          //
//                                                                      //
// HMPID class to perfom particle identification                        //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "AliHMPIDPid.h"       //class header
#include "AliHMPIDParam.h"     //class header
#include "AliHMPIDRecon.h"     //class header
#include <AliESDtrack.h>       //FindPid()
#include <TRandom.h>           //Resolution()

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliHMPIDPid::AliHMPIDPid():TNamed("HMPIDrec","HMPIDPid")
{
//..
//init of data members
//..
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Int_t nsp,Double_t *prob)
{
// Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
// from the given Cerenkov angle and momentum assuming no initial particle composition
// (i.e. apriory probability to be the particle of the given sort is the same for all sorts)

  AliPID *pPid = new AliPID();
  
  Double_t thetaCerExp = -999.;                                                                           
  if(pTrk->GetHMPIDsignal()<=0) thetaCerExp = pTrk->GetHMPIDsignal();                                                                           
  else                          thetaCerExp = pTrk->GetHMPIDsignal() - (Int_t)pTrk->GetHMPIDsignal();     //  measured thetaCherenkov
  
  if(thetaCerExp<=0){                                                                                     //HMPID does not find anything reasonable for this track, assign 0.2 for all species
    for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
    delete pPid ; pPid=0x0; return;
  } 
  
  Double_t p[3] = {0}, pmod = 0;
  if(pTrk->GetOuterHmpPxPyPz(p))  pmod = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);  // Momentum of the charged particle
  
  else {                                         
    for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
    delete pPid ; pPid=0x0; return;
  } 
  
  Double_t hTot=0;                               // Initialize the total height of the amplitude method
  Double_t *h = new Double_t [nsp];              // number of charged particles to be considered

  Bool_t desert = kTRUE;                                                                                                     //  Flag to evaluate if ThetaC is far ("desert") from the given Gaussians
  
  for(Int_t iPart=0;iPart<nsp;iPart++){                                                                                      //  for each particle
    
    h[iPart] = 0;                                                                                                            //  reset the height
    Double_t mass = pPid->ParticleMass(iPart);                                                                             //  with the given mass
    Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(AliHMPIDParam::Instance()->MeanIdxRad()*pmod);                   //  evaluate the theor. Theta Cherenkov
    if(cosThetaTh>1) continue;                                                                                               //  no light emitted, zero height
    Double_t thetaCerTh = TMath::ACos(cosThetaTh);                                                                           //  theoretical Theta Cherenkov
    Double_t sigmaRing = Resolution(thetaCerTh,pTrk);
    
    if(sigmaRing==0) {
      for(Int_t jPart=0;jPart<nsp;jPart++) prob[jPart]=1.0/(Float_t)nsp;
      delete pPid ; pPid=0x0; delete [] h; return;
    } 
    
    if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE;                                                                //   
    h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
    hTot    +=h[iPart]; //total height of all theoretical heights for normalization
    
  }//species loop

  for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
    
    if(!desert) prob[iPart]=h[iPart]/hTot;
    else prob[iPart]=1.0/(Float_t)nsp;            //all theoretical values are far away from experemental one
    
  }
  
  delete [] h;
  delete pPid ; pPid=0x0;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
{
  AliHMPIDParam *pParam = AliHMPIDParam::Instance();
      
  AliHMPIDRecon rec;
  Float_t xRa,yRa,thRa,phRa;
  pTrk->GetHMPIDtrk(xRa,yRa,thRa,phRa);
  rec.SetTrack(xRa,yRa,thRa,phRa);
  Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
  Int_t nPhots = (Int_t)(21.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);

  Double_t sigmatot = 0;
  Int_t nTrks = 20;
  for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
    Double_t invSigma = 0;
    Int_t nPhotsAcc = 0;
    for(Int_t j=0;j<nPhots;j++){
      Double_t phi = gRandom->Rndm()*TMath::TwoPi();
      TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
      if(!pParam->IsInside(pos.X(),pos.Y())) continue;
      if(pParam->IsInDead(pos.X(),pos.Y())) continue;
      Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
      if(sigma2!=0) {
        invSigma += 1./sigma2;
        nPhotsAcc++;
      }
    }      
    if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);  
  }
  return sigmatot/nTrks;
}
Commit	Line	Data
3278403b	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	// //
	18	// AliHMPIDPid //
	19	// //
	20	// HMPID class to perfom particle identification //
	21	// //
	22	//////////////////////////////////////////////////////////////////////////
	23
	24	#include "AliHMPIDPid.h" //class header
	25	#include "AliHMPIDParam.h" //class header
ac66a50d	26	#include "AliHMPIDRecon.h" //class header
9eb30f08	27	#include <AliESDtrack.h> //FindPid()
ac66a50d	28	#include <TRandom.h> //Resolution()
3278403b	29
3278403b	30	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
f21fc003	31	AliHMPIDPid::AliHMPIDPid():TNamed("HMPIDrec","HMPIDPid")
3278403b	32	{
	33	//..
	34	//init of data members
	35	//..
	36	}
	37	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
dac53a45	38	void AliHMPIDPid::FindPid(AliESDtrack pTrk,Int_t nsp,Double_t prob)
3278403b	39	{
ac66a50d	40	// Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
3278403b	41	// from the given Cerenkov angle and momentum assuming no initial particle composition
	42	// (i.e. apriory probability to be the particle of the given sort is the same for all sorts)
	43
ac66a50d	44	AliPID *pPid = new AliPID();
ac66a50d	45
6ba13ae0	46	Double_t thetaCerExp = -999.;
	47	if(pTrk->GetHMPIDsignal()<=0) thetaCerExp = pTrk->GetHMPIDsignal();
	48	else thetaCerExp = pTrk->GetHMPIDsignal() - (Int_t)pTrk->GetHMPIDsignal(); // measured thetaCherenkov
	49
	50	if(thetaCerExp<=0){ //HMPID does not find anything reasonable for this track, assign 0.2 for all species
dac53a45	51	for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
89e82ce1	52	delete pPid ; pPid=0x0; return;
9eb30f08	53	}
122305ea	54
	55	Double_t p[3] = {0}, pmod = 0;
	56	if(pTrk->GetOuterHmpPxPyPz(p)) pmod = TMath::Sqrt(p[0]p[0]+p[1]p[1]+p[2]*p[2]); // Momentum of the charged particle
	57
	58	else {
	59	for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
89e82ce1	60	delete pPid ; pPid=0x0; return;
122305ea	61	}
122305ea	62
ac66a50d	63	Double_t hTot=0; // Initialize the total height of the amplitude method
ac66a50d	64	Double_t *h = new Double_t [nsp]; // number of charged particles to be considered
9eb30f08	65
ac66a50d	66	Bool_t desert = kTRUE; // Flag to evaluate if ThetaC is far ("desert") from the given Gaussians
	67
	68	for(Int_t iPart=0;iPart<nsp;iPart++){ // for each particle
	69
	70	h[iPart] = 0; // reset the height
	71	Double_t mass = pPid->ParticleMass(iPart); // with the given mass
	72	Double_t cosThetaTh = TMath::Sqrt(massmass+pmodpmod)/(AliHMPIDParam::Instance()->MeanIdxRad()*pmod); // evaluate the theor. Theta Cherenkov
	73	if(cosThetaTh>1) continue; // no light emitted, zero height
	74	Double_t thetaCerTh = TMath::ACos(cosThetaTh); // theoretical Theta Cherenkov
	75	Double_t sigmaRing = Resolution(thetaCerTh,pTrk);
	76
	77	if(sigmaRing==0) {
	78	for(Int_t jPart=0;jPart<nsp;jPart++) prob[jPart]=1.0/(Float_t)nsp;
aa8db9d6	79	delete pPid ; pPid=0x0; delete [] h; return;
ac66a50d	80	}
	81
	82	if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE; //
	83	h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
9eb30f08	84	hTot +=h[iPart]; //total height of all theoretical heights for normalization
ac66a50d	85
9eb30f08	86	}//species loop
9eb30f08	87
dac53a45	88	for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
ac66a50d	89
ac66a50d	90	if(!desert) prob[iPart]=h[iPart]/hTot;
dac53a45	91	else prob[iPart]=1.0/(Float_t)nsp; //all theoretical values are far away from experemental one
ac66a50d	92
9eb30f08	93	}
ac66a50d	94
dac53a45	95	delete [] h;
ac66a50d	96	delete pPid ; pPid=0x0;
3278403b	97	}
3278403b	98	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ac66a50d	99	Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
	100	{
	101	AliHMPIDParam *pParam = AliHMPIDParam::Instance();
	102
	103	AliHMPIDRecon rec;
	104	Float_t xRa,yRa,thRa,phRa;
	105	pTrk->GetHMPIDtrk(xRa,yRa,thRa,phRa);
	106	rec.SetTrack(xRa,yRa,thRa,phRa);
	107	Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
	108	Int_t nPhots = (Int_t)(21.TMath::Sin(thetaCerTh)TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
	109
	110	Double_t sigmatot = 0;
	111	Int_t nTrks = 20;
	112	for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
	113	Double_t invSigma = 0;
	114	Int_t nPhotsAcc = 0;
	115	for(Int_t j=0;j<nPhots;j++){
	116	Double_t phi = gRandom->Rndm()*TMath::TwoPi();
	117	TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
	118	if(!pParam->IsInside(pos.X(),pos.Y())) continue;
	119	if(pParam->IsInDead(pos.X(),pos.Y())) continue;
	120	Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
	121	if(sigma2!=0) {
	122	invSigma += 1./sigma2;
	123	nPhotsAcc++;
	124	}
	125	}
	126	if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);
	127	}
	128	return sigmatot/nTrks;
	129	}