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

//-----------------------------------------------------------------//
//                                                                 //
//           Implementation of the TOF PID class                   //
//      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch         //
//                                                                 //
//-----------------------------------------------------------------//

#include "TMath.h"
#include "AliLog.h"

#include "AliESDtrack.h"
#include "AliESDEvent.h"

#include "AliTOFpidESD.h"

ClassImp(AliTOFpidESD)

//_________________________________________________________________________
AliTOFpidESD::AliTOFpidESD(): 
  fSigma(0),
  fRange(0),
  fPmax(0),         // zero at 0.5 GeV/c for pp
  fTime0(0)
{
}
//_________________________________________________________________________
AliTOFpidESD::AliTOFpidESD(Double_t *param):
  fSigma(param[0]),
  fRange(param[1]),
  fPmax(0),          // zero at 0.5 GeV/c for pp
  fTime0(0)
{
  //
  //  The main constructor
  //
  //

  //fPmax=TMath::Exp(-0.5*3*3)/fSigma; // ~3 sigma at 0.5 GeV/c for PbPb 
}

//_________________________________________________________________________
Double_t 
AliTOFpidESD::GetMismatchProbability(Double_t p, Double_t mass) const {
  //
  // Returns the probability of mismatching 
  // assuming 1/(p*beta)^2 scaling
  //
  const Double_t km=0.5;                   // "reference" momentum (GeV/c)

  Double_t ref2=km*km*km*km/(km*km + mass*mass);// "reference" (p*beta)^2
  Double_t p2beta2=p*p*p*p/(p*p + mass*mass);

  return fPmax*ref2/p2beta2;
}

//_________________________________________________________________________
Int_t AliTOFpidESD::MakePID(AliESDEvent *event, Double_t timeZero)
{
  //
  //  This function calculates the "detector response" PID probabilities
  //                Just for a bare hint... 

  AliDebug(1,Form("TOF PID Parameters: Sigma (ps)= %f, Range= %f",fSigma,fRange));
  AliDebug(1,"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \n");

  fTime0=timeZero;
  return MakePID(event);  
}

//_________________________________________________________________________
Int_t AliTOFpidESD::MakePID(AliESDEvent *event)
{
  //
  //  This function calculates the "detector response" PID probabilities
  //                Just for a bare hint... 

  Int_t ntrk=event->GetNumberOfTracks();
  AliESDtrack **tracks=new AliESDtrack*[ntrk];

  Int_t i;
  for (i=0; i<ntrk; i++) {
    AliESDtrack *t=event->GetTrack(i);
    tracks[i]=t;
  }

  for (i=0; i<ntrk; i++) {
    AliESDtrack *t=tracks[i];
    if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
    if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;

    Double_t time[AliPID::kSPECIES];
    if (!ExpectedSignals(t,time,AliPID::kSPECIES)) continue;
    Double_t sigma[AliPID::kSPECIES];
    if (!ExpectedSigmas(t,sigma,AliPID::kSPECIES)) continue;

    AliDebug(2,Form("Expected TOF signals [ps]: %f %f %f %f %f",
		    GetExpectedSignal(t,AliPID::kElectron),
		    GetExpectedSignal(t,AliPID::kMuon),
		    GetExpectedSignal(t,AliPID::kPion),
		    GetExpectedSignal(t,AliPID::kKaon),
		    GetExpectedSignal(t,AliPID::kProton)
		    ));

    AliDebug(2,Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
		    GetExpectedSigma(t,AliPID::kElectron),
		    GetExpectedSigma(t,AliPID::kMuon),
		    GetExpectedSigma(t,AliPID::kPion),
		    GetExpectedSigma(t,AliPID::kKaon),
		    GetExpectedSigma(t,AliPID::kProton)
		    ));

    AliDebug(2,Form("Expected TOF std deviations [number of expected sigmas]: %f %f %f %f %f",
		    GetNumberOfSigmas(t,AliPID::kElectron),
		    GetNumberOfSigmas(t,AliPID::kMuon),
		    GetNumberOfSigmas(t,AliPID::kPion),
		    GetNumberOfSigmas(t,AliPID::kKaon),
		    GetNumberOfSigmas(t,AliPID::kProton)
		    ));

    Double_t tof = t->GetTOFsignal() - fTime0;

    Double_t p[AliPID::kSPECIES];
    Bool_t mismatch = kTRUE, heavy = kTRUE;
    for (Int_t j=0; j<AliPID::kSPECIES; j++) {
      Double_t sig = sigma[j];
      if (TMath::Abs(tof-time[j]) > fRange*sig) {
	p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
      } else
        p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;

      // Check the mismatching
      Double_t mass = AliPID::ParticleMass(j);
      Double_t pm = GetMismatchProbability(t->GetP(),mass);
      if (p[j]>pm) mismatch = kFALSE;

      // Check for particles heavier than (AliPID::kSPECIES - 1)
      if (tof < (time[j] + fRange*sig)) heavy=kFALSE;

    }

    if (mismatch)
       for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;

    t->SetTOFpid(p);

    if (heavy) t->ResetStatus(AliESDtrack::kTOFpid);    

  }

  delete[] tracks;
  
  return 0;
}

//_________________________________________________________________________
Bool_t AliTOFpidESD::ExpectedSignals(const AliESDtrack *t,
                                     Double_t s[],Int_t n) const
{
  //
  // Return the expected PID signals for the involved particle species
  //

  if (n > AliPID::kSPECIESN)          return kFALSE;
  if ( !t->IsOn(AliESDtrack::kTIME) ) return kFALSE;

  Double_t time[AliPID::kSPECIESN];
  t->GetIntegratedTimes(time);
  for (Int_t i=0; i<n; i++) s[i]=time[i];
  return kTRUE;
}

//_________________________________________________________________________
Bool_t AliTOFpidESD::ExpectedSigmas(const AliESDtrack *t,
                                     Double_t s[],Int_t n) const
{
  //
  // Return the expected sigma of PID signals for the involved
  // particle species.
  // This approximate (but reasonable) formula takes into account the
  // relative momentum resolution.
  //

  Double_t time[AliPID::kSPECIESN];
  if ( !ExpectedSignals(t,time,n) ) return kFALSE;

  Double_t mom = t->GetP();
  Double_t dpp = 0.01;      //mean relative pt resolution;
  if (mom>0.5) dpp = 0.01*mom;
  for (Int_t i=0; i<n; i++) {
    Double_t mass = AliPID::ParticleMass(i);
    Double_t sigma = dpp*time[i]/(1.+ mom*mom/(mass*mass));
    s[i] = TMath::Sqrt(sigma*sigma + fSigma*fSigma);
  }
  return kTRUE;  
}

//_________________________________________________________________________
Bool_t AliTOFpidESD::NumberOfSigmas(const AliESDtrack *t,
                                    Double_t s[],Int_t n) const
{
  //
  // Returns the deviation of the actual PID signal from the expected
  // signal, in units of expected sigmas.
  //

  Double_t time[AliPID::kSPECIESN];
  if ( !ExpectedSignals(t,time,n) ) return kFALSE;

  if ( !ExpectedSigmas(t,s,n) ) return kFALSE;
  
  Double_t tof = t->GetTOFsignal() - fTime0;
  for (Int_t i=0; i<n; i++) s[i] = (time[i]-tof)/s[i];
      
  return kTRUE;
}

//_________________________________________________________________________
 Double_t AliTOFpidESD::GetExpectedSignal(const AliESDtrack *t,
                                          AliPID::EParticleType n) const
{
  //
  // Return the expected PID signal for the specified particle type.
  // If the operation is not possible, return a negative value.
  //

  if (Int_t(n) >= AliPID::kSPECIESN)        return -1.;
  if ( !t->IsOn(AliESDtrack::kTIME) ) return -1.;

  Double_t time[AliPID::kSPECIESN];
  t->GetIntegratedTimes(time);

  return time[n];
}

//_________________________________________________________________________
 Double_t AliTOFpidESD::GetExpectedSigma(const AliESDtrack *t,
                                         AliPID::EParticleType n) const
{
  //
  // Return the expected sigma of the PID signal for the specified
  // particle type.
  // If the operation is not possible, return a negative value.
  //

  Double_t time[AliPID::kSPECIESN];
  if ( !ExpectedSignals(t,time,AliPID::kSPECIESN) ) return -1.;

  Double_t mom = t->GetP();
  Double_t dpp = 0.01;      //mean relative pt resolution;
  if (mom>0.5) dpp = 0.01*mom;

  Double_t mass = AliPID::ParticleMass(n);
  Double_t sigma = dpp*time[n]/(1.+ mom*mom/(mass*mass));

  return TMath::Sqrt(sigma*sigma + fSigma*fSigma);
}

//_________________________________________________________________________
 Double_t AliTOFpidESD::GetNumberOfSigmas(const AliESDtrack *t,
                                          AliPID::EParticleType n) const
{
  //
  // Returns the deviation of the actual PID signal from the expected
  // signal for the specified particle type, in units of expected
  // sigmas.
  // If the operation is not possible, return a negative value.
  //

  Double_t time=GetExpectedSignal(t,n);;
  if (time < 0.) return -1.;

  Double_t sigma=GetExpectedSigma(t,n);
  if (sigma < 0.) return -1;
  
  Double_t tof=t->GetTOFsignal() - fTime0;
  return (time-tof)/sigma;
}
Commit	Line	Data
c630aafd	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
0e46b9ae	16	//-----------------------------------------------------------------//
	17	// //
	18	// Implementation of the TOF PID class //
	19	// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
	20	// //
	21	//-----------------------------------------------------------------//
	22
	23	#include "TMath.h"
02cdc6ac	24	#include "AliLog.h"
0e46b9ae	25
c630aafd	26	#include "AliESDtrack.h"
af885e0f	27	#include "AliESDEvent.h"
c630aafd	28
0e46b9ae	29	#include "AliTOFpidESD.h"
0e46b9ae	30
c630aafd	31	ClassImp(AliTOFpidESD)
c630aafd	32
0f4a7374	33	//_________________________________________________________________________
7a8614f3	34	AliTOFpidESD::AliTOFpidESD():
655e379f	35	fSigma(0),
7a8614f3	36	fRange(0),
cc79c586	37	fPmax(0), // zero at 0.5 GeV/c for pp
cc79c586	38	fTime0(0)
655e379f	39	{
	40	}
	41	//_________________________________________________________________________
	42	AliTOFpidESD::AliTOFpidESD(Double_t *param):
7a8614f3	43	fSigma(param[0]),
7a8614f3	44	fRange(param[1]),
cc79c586	45	fPmax(0), // zero at 0.5 GeV/c for pp
cc79c586	46	fTime0(0)
7a8614f3	47	{
0f4a7374	48	//
	49	// The main constructor
	50	//
7a8614f3	51	//
	52
	53	//fPmax=TMath::Exp(-0.533)/fSigma; // ~3 sigma at 0.5 GeV/c for PbPb
	54	}
	55
cc79c586	56	//_________________________________________________________________________
7a8614f3	57	Double_t
	58	AliTOFpidESD::GetMismatchProbability(Double_t p, Double_t mass) const {
	59	//
	60	// Returns the probability of mismatching
	61	// assuming 1/(p*beta)^2 scaling
	62	//
62bdc654	63	const Double_t km=0.5; // "reference" momentum (GeV/c)
7a8614f3	64
62bdc654	65	Double_t ref2=kmkmkmkm/(kmkm + massmass);// "reference" (pbeta)^2
7a8614f3	66	Double_t p2beta2=pppp/(pp + mass*mass);
0f4a7374	67
7a8614f3	68	return fPmax*ref2/p2beta2;
0f4a7374	69	}
0f4a7374	70
bc9f08da	71	//_________________________________________________________________________
af885e0f	72	Int_t AliTOFpidESD::MakePID(AliESDEvent *event, Double_t timeZero)
bc9f08da	73	{
	74	//
	75	// This function calculates the "detector response" PID probabilities
	76	// Just for a bare hint...
	77
02cdc6ac	78	AliDebug(1,Form("TOF PID Parameters: Sigma (ps)= %f, Range= %f",fSigma,fRange));
02cdc6ac	79	AliDebug(1,"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \n");
7a8614f3	80
cc79c586	81	fTime0=timeZero;
	82	return MakePID(event);
	83	}
	84
	85	//_________________________________________________________________________
	86	Int_t AliTOFpidESD::MakePID(AliESDEvent *event)
	87	{
	88	//
	89	// This function calculates the "detector response" PID probabilities
	90	// Just for a bare hint...
	91
bc9f08da	92	Int_t ntrk=event->GetNumberOfTracks();
	93	AliESDtrack *tracks=new AliESDtrack[ntrk];
	94
	95	Int_t i;
	96	for (i=0; i<ntrk; i++) {
	97	AliESDtrack *t=event->GetTrack(i);
	98	tracks[i]=t;
	99	}
	100
	101	for (i=0; i<ntrk; i++) {
	102	AliESDtrack *t=tracks[i];
	103	if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
	104	if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
cc79c586	105
	106	Double_t time[AliPID::kSPECIES];
	107	if (!ExpectedSignals(t,time,AliPID::kSPECIES)) continue;
	108	Double_t sigma[AliPID::kSPECIES];
	109	if (!ExpectedSigmas(t,sigma,AliPID::kSPECIES)) continue;
	110
	111	AliDebug(2,Form("Expected TOF signals [ps]: %f %f %f %f %f",
	112	GetExpectedSignal(t,AliPID::kElectron),
	113	GetExpectedSignal(t,AliPID::kMuon),
	114	GetExpectedSignal(t,AliPID::kPion),
	115	GetExpectedSignal(t,AliPID::kKaon),
	116	GetExpectedSignal(t,AliPID::kProton)
	117	));
	118
	119	AliDebug(2,Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
	120	GetExpectedSigma(t,AliPID::kElectron),
	121	GetExpectedSigma(t,AliPID::kMuon),
	122	GetExpectedSigma(t,AliPID::kPion),
	123	GetExpectedSigma(t,AliPID::kKaon),
	124	GetExpectedSigma(t,AliPID::kProton)
	125	));
	126
	127	AliDebug(2,Form("Expected TOF std deviations [number of expected sigmas]: %f %f %f %f %f",
	128	GetNumberOfSigmas(t,AliPID::kElectron),
	129	GetNumberOfSigmas(t,AliPID::kMuon),
	130	GetNumberOfSigmas(t,AliPID::kPion),
	131	GetNumberOfSigmas(t,AliPID::kKaon),
	132	GetNumberOfSigmas(t,AliPID::kProton)
	133	));
	134
	135	Double_t tof = t->GetTOFsignal() - fTime0;
	136
	137	Double_t p[AliPID::kSPECIES];
	138	Bool_t mismatch = kTRUE, heavy = kTRUE;
bc9f08da	139	for (Int_t j=0; j<AliPID::kSPECIES; j++) {
cc79c586	140	Double_t sig = sigma[j];
	141	if (TMath::Abs(tof-time[j]) > fRange*sig) {
	142	p[j] = TMath::Exp(-0.5fRangefRange)/sig;
	143	} else
	144	p[j] = TMath::Exp(-0.5(tof-time[j])(tof-time[j])/(sig*sig))/sig;
7a8614f3	145
7a8614f3	146	// Check the mismatching
cc79c586	147	Double_t mass = AliPID::ParticleMass(j);
	148	Double_t pm = GetMismatchProbability(t->GetP(),mass);
	149	if (p[j]>pm) mismatch = kFALSE;
7a8614f3	150
7a8614f3	151	// Check for particles heavier than (AliPID::kSPECIES - 1)
cc79c586	152	if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
7a8614f3	153
bc9f08da	154	}
7a8614f3	155
	156	if (mismatch)
	157	for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
	158
bc9f08da	159	t->SetTOFpid(p);
7a8614f3	160
	161	if (heavy) t->ResetStatus(AliESDtrack::kTOFpid);
	162
bc9f08da	163	}
	164
	165	delete[] tracks;
	166
	167	return 0;
	168	}
	169
c630aafd	170	//_________________________________________________________________________
cc79c586	171	Bool_t AliTOFpidESD::ExpectedSignals(const AliESDtrack *t,
cc79c586	172	Double_t s[],Int_t n) const
c630aafd	173	{
c630aafd	174	//
cc79c586	175	// Return the expected PID signals for the involved particle species
cc79c586	176	//
c630aafd	177
cc79c586	178	if (n > AliPID::kSPECIESN) return kFALSE;
cc79c586	179	if ( !t->IsOn(AliESDtrack::kTIME) ) return kFALSE;
3f83f224	180
cc79c586	181	Double_t time[AliPID::kSPECIESN];
	182	t->GetIntegratedTimes(time);
	183	for (Int_t i=0; i<n; i++) s[i]=time[i];
	184	return kTRUE;
	185	}
	186
	187	//_________________________________________________________________________
	188	Bool_t AliTOFpidESD::ExpectedSigmas(const AliESDtrack *t,
	189	Double_t s[],Int_t n) const
	190	{
	191	//
	192	// Return the expected sigma of PID signals for the involved
	193	// particle species.
	194	// This approximate (but reasonable) formula takes into account the
	195	// relative momentum resolution.
	196	//
	197
	198	Double_t time[AliPID::kSPECIESN];
	199	if ( !ExpectedSignals(t,time,n) ) return kFALSE;
	200
	201	Double_t mom = t->GetP();
	202	Double_t dpp = 0.01; //mean relative pt resolution;
	203	if (mom>0.5) dpp = 0.01*mom;
	204	for (Int_t i=0; i<n; i++) {
	205	Double_t mass = AliPID::ParticleMass(i);
	206	Double_t sigma = dpptime[i]/(1.+ mommom/(mass*mass));
	207	s[i] = TMath::Sqrt(sigmasigma + fSigmafSigma);
3f83f224	208	}
cc79c586	209	return kTRUE;
cc79c586	210	}
3f83f224	211
cc79c586	212	//_________________________________________________________________________
	213	Bool_t AliTOFpidESD::NumberOfSigmas(const AliESDtrack *t,
	214	Double_t s[],Int_t n) const
	215	{
	216	//
	217	// Returns the deviation of the actual PID signal from the expected
	218	// signal, in units of expected sigmas.
	219	//
7a8614f3	220
cc79c586	221	Double_t time[AliPID::kSPECIESN];
cc79c586	222	if ( !ExpectedSignals(t,time,n) ) return kFALSE;
7a8614f3	223
cc79c586	224	if ( !ExpectedSigmas(t,s,n) ) return kFALSE;
	225
	226	Double_t tof = t->GetTOFsignal() - fTime0;
	227	for (Int_t i=0; i<n; i++) s[i] = (time[i]-tof)/s[i];
	228
	229	return kTRUE;
	230	}
7a8614f3	231
cc79c586	232	//_________________________________________________________________________
	233	Double_t AliTOFpidESD::GetExpectedSignal(const AliESDtrack *t,
	234	AliPID::EParticleType n) const
	235	{
	236	//
	237	// Return the expected PID signal for the specified particle type.
	238	// If the operation is not possible, return a negative value.
	239	//
7a8614f3	240
cc79c586	241	if (Int_t(n) >= AliPID::kSPECIESN) return -1.;
cc79c586	242	if ( !t->IsOn(AliESDtrack::kTIME) ) return -1.;
7a8614f3	243
cc79c586	244	Double_t time[AliPID::kSPECIESN];
cc79c586	245	t->GetIntegratedTimes(time);
7a8614f3	246
cc79c586	247	return time[n];
cc79c586	248	}
7a8614f3	249
cc79c586	250	//_________________________________________________________________________
	251	Double_t AliTOFpidESD::GetExpectedSigma(const AliESDtrack *t,
	252	AliPID::EParticleType n) const
	253	{
	254	//
	255	// Return the expected sigma of the PID signal for the specified
	256	// particle type.
	257	// If the operation is not possible, return a negative value.
	258	//
c630aafd	259
cc79c586	260	Double_t time[AliPID::kSPECIESN];
	261	if ( !ExpectedSignals(t,time,AliPID::kSPECIESN) ) return -1.;
	262
	263	Double_t mom = t->GetP();
	264	Double_t dpp = 0.01; //mean relative pt resolution;
	265	if (mom>0.5) dpp = 0.01*mom;
	266
	267	Double_t mass = AliPID::ParticleMass(n);
	268	Double_t sigma = dpptime[n]/(1.+ mommom/(mass*mass));
	269
	270	return TMath::Sqrt(sigmasigma + fSigmafSigma);
c630aafd	271	}
0f4a7374	272
cc79c586	273	//_________________________________________________________________________
	274	Double_t AliTOFpidESD::GetNumberOfSigmas(const AliESDtrack *t,
	275	AliPID::EParticleType n) const
	276	{
	277	//
	278	// Returns the deviation of the actual PID signal from the expected
	279	// signal for the specified particle type, in units of expected
	280	// sigmas.
	281	// If the operation is not possible, return a negative value.
	282	//
	283
	284	Double_t time=GetExpectedSignal(t,n);;
	285	if (time < 0.) return -1.;
	286
	287	Double_t sigma=GetExpectedSigma(t,n);
	288	if (sigma < 0.) return -1;
	289
	290	Double_t tof=t->GetTOFsignal() - fTime0;
	291	return (time-tof)/sigma;
	292	}