[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
// Very naive one... Should be made better by the detector experts...
//      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
//-----------------------------------------------------------------
#include "TFile.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TError.h"

#include "AliTOFpidESD.h"
#include "AliESD.h"
#include "AliESDtrack.h"
#include "AliTOFdigit.h"
#include "AliTOFGeometry.h"

#include <stdlib.h>


ClassImp(AliTOFpidESD)

//_________________________________________________________________________
AliTOFpidESD::AliTOFpidESD(Double_t *param) {
  //
  //  The main constructor
  //
  fR=378.; 
  fDy=AliTOFGeometry::XPad(); fDz=AliTOFGeometry::ZPad(); 
  fN=0; fEventN=0;

  fSigma=param[0];
  fRange=param[1];

}

//_________________________________________________________________________
Int_t AliTOFpidESD::LoadClusters(TTree *dTree, AliTOFGeometry *geom) {
  //--------------------------------------------------------------------
  //This function loads the TOF clusters
  //--------------------------------------------------------------------
  TBranch *branch=dTree->GetBranch("TOF");
  if (!branch) { 
    Error("LoadClusters"," can't get the branch with the TOF digits !\n");
    return 1;
  }

  TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
  branch->SetAddress(&digits);

  dTree->GetEvent(0);
  Int_t nd=digits->GetEntriesFast();
  Info("LoadClusters","number of digits: %d",nd);

  for (Int_t i=0; i<nd; i++) {
    AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
    Int_t dig[5]; Float_t g[3];
    dig[0]=d->GetSector();
    dig[1]=d->GetPlate();
    dig[2]=d->GetStrip();
    dig[3]=d->GetPadz();
    dig[4]=d->GetPadx();

    geom->GetPos(dig,g);

    Double_t h[5];
    h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
    h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; 
    h[3]=d->GetTdc(); h[4]=d->GetAdc();

    AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
    InsertCluster(cl);
  }  

  return 0;
}

//_________________________________________________________________________
void AliTOFpidESD::UnloadClusters() {
  //--------------------------------------------------------------------
  //This function unloads TOF clusters
  //--------------------------------------------------------------------
  for (Int_t i=0; i<fN; i++) delete fClusters[i];
  fN=0;
}

//_________________________________________________________________________
Int_t AliTOFpidESD::InsertCluster(AliTOFcluster *c) {
  //--------------------------------------------------------------------
  //This function adds a cluster to the array of clusters sorted in Z
  //--------------------------------------------------------------------
  if (fN==kMaxCluster) {
    Error("InsertCluster","Too many clusters !\n");
    return 1;
  }

  if (fN==0) {fClusters[fN++]=c; return 0;}
  Int_t i=FindClusterIndex(c->GetZ());
  memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTOFcluster*));
  fClusters[i]=c; fN++;

  return 0;
}

//_________________________________________________________________________
Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
  //--------------------------------------------------------------------
  // This function returns the index of the nearest cluster 
  //--------------------------------------------------------------------
  if (fN==0) return 0;
  if (z <= fClusters[0]->GetZ()) return 0;
  if (z > fClusters[fN-1]->GetZ()) return fN;
  Int_t b=0, e=fN-1, m=(b+e)/2;
  for (; b<e; m=(b+e)/2) {
    if (z > fClusters[m]->GetZ()) b=m+1;
    else e=m; 
  }
  return m;
}

static int cmp(const void *p1, const void *p2) {
  AliESDtrack *t1=*((AliESDtrack**)p1);
  AliESDtrack *t2=*((AliESDtrack**)p2);
  Double_t c1[15]; t1->GetExternalCovariance(c1);
  Double_t c2[15]; t2->GetExternalCovariance(c2);
  if (c1[0]*c1[2] <c2[0]*c2[2]) return -1;
  if (c1[0]*c1[2]==c2[0]*c2[2]) return 0;
  return 1;
}

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

  static const Double_t kMasses[]={
    0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
  };

  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;
  }
  qsort(tracks,ntrk,sizeof(AliESDtrack*),cmp);

  Int_t nmatch=0;
  for (i=0; i<ntrk; i++) {
    AliESDtrack *t=tracks[i];

    if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
    if ((t->GetStatus()&AliESDtrack::kTRDStop)!=0) continue;

    Double_t time[10]; t->GetIntegratedTimes(time);

    Double_t x,par[5]; t->GetExternalParameters(x,par);
    Double_t cov[15]; t->GetExternalCovariance(cov);

Double_t dphi=(5*TMath::Sqrt(cov[0]) + 0.5*fDy + 2.5*TMath::Abs(par[2]))/fR; 
Double_t dz=5*TMath::Sqrt(cov[2]) + 0.5*fDz + 2.5*TMath::Abs(par[3]);

    Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha();
    if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
    if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
    Double_t z=par[1];

    Double_t d2max=1000.;
    Int_t index=-1;
    for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
      AliTOFcluster *c=fClusters[k];

      if (c->GetZ() > z+dz) break;
      if (c->IsUsed()) continue;

      Double_t tof=50*c->GetTDC()+32;
      if (t->GetIntegratedLength()/tof > 0.031) continue;
      if (tof>35000) continue;

      Double_t dph=TMath::Abs(c->GetPhi()-phi);
      if (dph>TMath::Pi()) dph=2*TMath::Pi()-dph; //Thanks to B.Zagreev
      if (dph>dphi) continue;

      Double_t d2=dph*dph*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
      if (d2 > d2max) continue;

      d2max=d2;
      index=k;
    }

    if (index<0) {
      //Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel()));
       continue;
    }

    nmatch++;

    AliTOFcluster *c=fClusters[index];
    c->Use();

    Double_t tof=50*c->GetTDC()+32; // in ps
    t->SetTOFsignal(tof);
    t->SetTOFcluster(c->GetIndex());

    if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;

    //track length correction
    Double_t rc=TMath::Sqrt(c->GetR()*c->GetR() + c->GetZ()*c->GetZ());
    Double_t rt=TMath::Sqrt(x*x + par[0]*par[0] + par[1]*par[1]);
    Double_t dlt=rc-rt;

    Double_t p[10];
    Double_t mom=t->GetP();
    for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
      Double_t mass=kMasses[j];
      Double_t dpp=0.01;      //mean relative pt resolution;
      if (mom>0.5) dpp=0.01*mom;
      Double_t sigma=dpp*time[j]/(1.+ mom*mom/(mass*mass));
      sigma=TMath::Sqrt(sigma*sigma + fSigma*fSigma);

      time[j]+=dlt/3e-2*TMath::Sqrt(mom*mom+mass*mass)/mom;

      if (TMath::Abs(tof-time[j]) > fRange*sigma) {
	p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
        continue;
      }
      p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma))/sigma;
    }
    /*
    if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
       cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
       continue;
    }
    */
    t->SetTOFpid(p);

  }

  Info("MakePID","Number of matched ESD track: %d",nmatch);

  delete[] tracks;

  return 0;
}
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
	16	//-----------------------------------------------------------------
	17	// Implementation of the TOF PID class
	18	// Very naive one... Should be made better by the detector experts...
	19	// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
	20	//-----------------------------------------------------------------
	21	#include "TFile.h"
	22	#include "TTree.h"
	23	#include "TClonesArray.h"
3f83f224	24	#include "TError.h"
c630aafd	25
	26	#include "AliTOFpidESD.h"
	27	#include "AliESD.h"
	28	#include "AliESDtrack.h"
	29	#include "AliTOFdigit.h"
8c29135e	30	#include "AliTOFGeometry.h"
c630aafd	31
3f83f224	32	#include <stdlib.h>
3f83f224	33
c630aafd	34
	35	ClassImp(AliTOFpidESD)
	36
0f4a7374	37	//_________________________________________________________________________
8c29135e	38	AliTOFpidESD::AliTOFpidESD(Double_t *param) {
0f4a7374	39	//
	40	// The main constructor
	41	//
0f4a7374	42	fR=378.;
8c29135e	43	fDy=AliTOFGeometry::XPad(); fDz=AliTOFGeometry::ZPad();
0f4a7374	44	fN=0; fEventN=0;
	45
	46	fSigma=param[0];
	47	fRange=param[1];
	48
	49	}
	50
0f4a7374	51	//_________________________________________________________________________
8c29135e	52	Int_t AliTOFpidESD::LoadClusters(TTree dTree, AliTOFGeometry geom) {
c630aafd	53	//--------------------------------------------------------------------
	54	//This function loads the TOF clusters
	55	//--------------------------------------------------------------------
	56	TBranch *branch=dTree->GetBranch("TOF");
	57	if (!branch) {
	58	Error("LoadClusters"," can't get the branch with the TOF digits !\n");
	59	return 1;
	60	}
	61
	62	TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
	63	branch->SetAddress(&digits);
	64
	65	dTree->GetEvent(0);
	66	Int_t nd=digits->GetEntriesFast();
	67	Info("LoadClusters","number of digits: %d",nd);
	68
	69	for (Int_t i=0; i<nd; i++) {
	70	AliTOFdigit d=(AliTOFdigit)digits->UncheckedAt(i);
	71	Int_t dig[5]; Float_t g[3];
	72	dig[0]=d->GetSector();
	73	dig[1]=d->GetPlate();
	74	dig[2]=d->GetStrip();
	75	dig[3]=d->GetPadz();
	76	dig[4]=d->GetPadx();
	77
8c29135e	78	geom->GetPos(dig,g);
c630aafd	79
	80	Double_t h[5];
	81	h[0]=TMath::Sqrt(g[0]g[0]+g[1]g[1]);
	82	h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
	83	h[3]=d->GetTdc(); h[4]=d->GetAdc();
	84
3f83f224	85	AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
c630aafd	86	InsertCluster(cl);
	87	}
	88
	89	return 0;
	90	}
	91
0f4a7374	92	//_________________________________________________________________________
c630aafd	93	void AliTOFpidESD::UnloadClusters() {
	94	//--------------------------------------------------------------------
	95	//This function unloads TOF clusters
	96	//--------------------------------------------------------------------
	97	for (Int_t i=0; i<fN; i++) delete fClusters[i];
	98	fN=0;
	99	}
	100
0f4a7374	101	//_________________________________________________________________________
c630aafd	102	Int_t AliTOFpidESD::InsertCluster(AliTOFcluster *c) {
	103	//--------------------------------------------------------------------
	104	//This function adds a cluster to the array of clusters sorted in Z
	105	//--------------------------------------------------------------------
	106	if (fN==kMaxCluster) {
	107	Error("InsertCluster","Too many clusters !\n");
	108	return 1;
	109	}
	110
	111	if (fN==0) {fClusters[fN++]=c; return 0;}
	112	Int_t i=FindClusterIndex(c->GetZ());
	113	memmove(fClusters+i+1 ,fClusters+i,(fN-i)sizeof(AliTOFcluster));
	114	fClusters[i]=c; fN++;
	115
	116	return 0;
	117	}
	118
0f4a7374	119	//_________________________________________________________________________
c630aafd	120	Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
	121	//--------------------------------------------------------------------
	122	// This function returns the index of the nearest cluster
	123	//--------------------------------------------------------------------
	124	if (fN==0) return 0;
	125	if (z <= fClusters[0]->GetZ()) return 0;
	126	if (z > fClusters[fN-1]->GetZ()) return fN;
	127	Int_t b=0, e=fN-1, m=(b+e)/2;
	128	for (; b<e; m=(b+e)/2) {
	129	if (z > fClusters[m]->GetZ()) b=m+1;
	130	else e=m;
	131	}
	132	return m;
	133	}
	134
3f83f224	135	static int cmp(const void p1, const void p2) {
	136	AliESDtrack t1=((AliESDtrack**)p1);
	137	AliESDtrack t2=((AliESDtrack**)p2);
	138	Double_t c1[15]; t1->GetExternalCovariance(c1);
	139	Double_t c2[15]; t2->GetExternalCovariance(c2);
	140	if (c1[0]c1[2] <c2[0]c2[2]) return -1;
	141	if (c1[0]c1[2]==c2[0]c2[2]) return 0;
	142	return 1;
	143	}
	144
c630aafd	145	//_________________________________________________________________________
	146	Int_t AliTOFpidESD::MakePID(AliESD *event)
	147	{
	148	//
	149	// This function calculates the "detector response" PID probabilities
	150	// Just for a bare hint...
	151
3f83f224	152	static const Double_t kMasses[]={
c630aafd	153	0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
	154	};
	155
	156	Int_t ntrk=event->GetNumberOfTracks();
3f83f224	157	AliESDtrack *tracks=new AliESDtrack[ntrk];
	158
	159	Int_t i;
	160	for (i=0; i<ntrk; i++) {
c630aafd	161	AliESDtrack *t=event->GetTrack(i);
3f83f224	162	tracks[i]=t;
	163	}
	164	qsort(tracks,ntrk,sizeof(AliESDtrack*),cmp);
	165
	166	Int_t nmatch=0;
	167	for (i=0; i<ntrk; i++) {
	168	AliESDtrack *t=tracks[i];
c630aafd	169
c630aafd	170	if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
b861b2b9	171	if ((t->GetStatus()&AliESDtrack::kTRDStop)!=0) continue;
c630aafd	172
c5957288	173	Double_t time[10]; t->GetIntegratedTimes(time);
c5957288	174
c630aafd	175	Double_t x,par[5]; t->GetExternalParameters(x,par);
	176	Double_t cov[15]; t->GetExternalCovariance(cov);
	177
3f83f224	178	Double_t dphi=(5TMath::Sqrt(cov[0]) + 0.5fDy + 2.5*TMath::Abs(par[2]))/fR;
3f83f224	179	Double_t dz=5TMath::Sqrt(cov[2]) + 0.5fDz + 2.5*TMath::Abs(par[3]);
c630aafd	180
	181	Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha();
	182	if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
	183	if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
	184	Double_t z=par[1];
	185
	186	Double_t d2max=1000.;
	187	Int_t index=-1;
	188	for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
	189	AliTOFcluster *c=fClusters[k];
	190
	191	if (c->GetZ() > z+dz) break;
	192	if (c->IsUsed()) continue;
	193
c5957288	194	Double_t tof=50*c->GetTDC()+32;
	195	if (t->GetIntegratedLength()/tof > 0.031) continue;
	196	if (tof>35000) continue;
	197
3f83f224	198	Double_t dph=TMath::Abs(c->GetPhi()-phi);
c5957288	199	if (dph>TMath::Pi()) dph=2*TMath::Pi()-dph; //Thanks to B.Zagreev
3f83f224	200	if (dph>dphi) continue;
c630aafd	201
3f83f224	202	Double_t d2=dphdphfRfR + (c->GetZ()-z)(c->GetZ()-z);
c630aafd	203	if (d2 > d2max) continue;
	204
	205	d2max=d2;
	206	index=k;
	207	}
	208
	209	if (index<0) {
	210	//Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel()));
	211	continue;
	212	}
	213
	214	nmatch++;
	215
	216	AliTOFcluster *c=fClusters[index];
3f83f224	217	c->Use();
	218
	219	Double_t tof=50*c->GetTDC()+32; // in ps
	220	t->SetTOFsignal(tof);
	221	t->SetTOFcluster(c->GetIndex());
	222
	223	if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
c630aafd	224
3f83f224	225	//track length correction
	226	Double_t rc=TMath::Sqrt(c->GetR()c->GetR() + c->GetZ()c->GetZ());
	227	Double_t rt=TMath::Sqrt(xx + par[0]par[0] + par[1]*par[1]);
	228	Double_t dlt=rc-rt;
c630aafd	229
	230	Double_t p[10];
	231	Double_t mom=t->GetP();
	232	for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
3f83f224	233	Double_t mass=kMasses[j];
	234	Double_t dpp=0.01; //mean relative pt resolution;
	235	if (mom>0.5) dpp=0.01*mom;
	236	Double_t sigma=dpptime[j]/(1.+ mommom/(mass*mass));
c630aafd	237	sigma=TMath::Sqrt(sigmasigma + fSigmafSigma);
3f83f224	238
	239	time[j]+=dlt/3e-2TMath::Sqrt(mommom+mass*mass)/mom;
	240
c630aafd	241	if (TMath::Abs(tof-time[j]) > fRange*sigma) {
431be10d	242	p[j]=TMath::Exp(-0.5fRangefRange)/sigma;
c630aafd	243	continue;
c630aafd	244	}
431be10d	245	p[j]=TMath::Exp(-0.5(tof-time[j])(tof-time[j])/(sigma*sigma))/sigma;
c630aafd	246	}
3f83f224	247	/*
	248	if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
	249	cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
	250	continue;
	251	}
	252	*/
c630aafd	253	t->SetTOFpid(p);
c630aafd	254
c630aafd	255	}
	256
	257	Info("MakePID","Number of matched ESD track: %d",nmatch);
	258
3f83f224	259	delete[] tracks;
3f83f224	260
c630aafd	261	return 0;
c630aafd	262	}
0f4a7374	263