[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 <TGeant3.h>              // For DtoM and InitGeo, should  
#include "../STRUCT/AliBODY.h"    // be removed after we switch to 
#include "../STRUCT/AliFRAMEv2.h" // AliTOGeometry::GetPos
#include "AliTOFv2FHoles.h"

#include <stdlib.h>


ClassImp(AliTOFpidESD)

//_________________________________________________________________________
static Int_t InitGeo() {
  //gSystem->Load("libgeant321");
  //new TGeant3("C++ Interface to Geant3");

  AliBODY *body = new AliBODY("BODY", "Alice envelop");
  body->CreateGeometry();

  AliFRAMEv2 *frame = new AliFRAMEv2("FRAME", "Space Frame");
  frame->SetHoles(1);
  frame->CreateGeometry();

  AliTOF *tof = new AliTOFv2FHoles("TOF", "TOF with Holes");
  tof->CreateGeometry();

  return 0;
} 

//_________________________________________________________________________
static Int_t DtoM(Int_t *dig, Float_t *g) {
  const Int_t kMAX=13;
  Int_t lnam[kMAX],lnum[kMAX];

  TGeant3 *geant3=(TGeant3*)gMC;
  if (!geant3) {::Info("DtoM","no geant3 found !"); return 1;}

  strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1;
  strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1;

//4 padx  if z<=0 then ...
if ((dig[1]==4)||(dig[1]==3)) dig[4]=AliTOFGeometry::NpadX()-dig[4];
else if (dig[1]==2) {
   if (dig[2]>7) dig[4]=AliTOFGeometry::NpadX()-dig[4];
   else if (dig[2]==7) {
      if (dig[3]==1) dig[4]=AliTOFGeometry::NpadX()-dig[4];
      else dig[4]+=1; 
   } else dig[4]+=1; 
} else dig[4]+=1;

//3 padz
if ((dig[1]==3)||(dig[1]==4)) dig[3]=AliTOFGeometry::NpadZ()-dig[3];
else dig[3]+=1;

//2 strip
if (dig[1]==0) dig[2]=AliTOFGeometry::NStripC()-dig[2]; 
else if (dig[1]==1) dig[2]=AliTOFGeometry::NStripB()-dig[2];
else dig[2]+=1; 

//1 module
dig[1]=AliTOFGeometry::NPlates()-dig[1];

//0 sector
dig[0]+=(AliTOFGeometry::NSectors()/2); dig[0]%=AliTOFGeometry::NSectors();
dig[0]+=1;

  //sector
  switch (dig[0]) {
     case 1:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=9;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 2:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=10;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 3:
         strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=4;
         strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
         if (dig[1]>3) lnum[3]=2;
         break;
     case 4:
         strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=5;
         strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
         if (dig[1]>3) lnum[3]=2;
         break;
     case 5:
         strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=1;
         strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
         if (dig[1]>3) lnum[3]=2;
         break;
     case 6:
         strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=2;
         strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
         if (dig[1]>3) lnum[3]=2;
         break;
     case 7:
         strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=3;
         strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
         if (dig[1]>3) lnum[3]=2;
         break;
     case 8:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=1;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 9:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=2;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 10:
         strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=1;
         strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
         if (dig[1]>4) lnum[3]=2;
         break;
     case 11:
         strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=2;
         strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
         if (dig[1]>4) lnum[3]=2;
         break;
     case 12:
         strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=3;
         strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
         if (dig[1]>4) lnum[3]=2;
         break;
     case 13:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=3;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 14:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=4;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 15:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=5;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 16:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=6;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 17:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=7;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     case 18:
         strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=8;
         strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
         break;
     default:
         ::Info("DtoM","Wrong sector number : %d !",dig[0]);
         return 2;
  }
  //module
  switch (dig[1]) {
     case 1:
        strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
        if (dig[0]==1  || dig[0]==2 ) lnum[4]=2;
        if (dig[0]>=8  && dig[0]<=9) lnum[4]=2;
        if (dig[0]>=13 && dig[0]<=18) lnum[4]=2;
        strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
        break;
     case 2:
        strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
        if (dig[0]==1  || dig[0]==2 ) lnum[4]=2;
        if (dig[0]>=8  && dig[0]<=9) lnum[4]=2;
        if (dig[0]>=13 && dig[0]<=18) lnum[4]=2;
        strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
        break;
     case 3:
        strncpy((Char_t*)(lnam+4),"FTOA",4); lnum[4]=0;
        strncpy((Char_t*)(lnam+5),"FLTA",4); lnum[5]=0;
        break;
     case 4:
        strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
        strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
        break;
     case 5:
        strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
        strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
        break;
     default:
        ::Info("DtoM","Wrong module number : %d !",dig[1]);
        return 3;
  }

  //strip
  if (dig[2]<0/*1*/ || dig[2]>19) {
     ::Info("DtoM","Wrong strip number : %d !",dig[2]);
     return 3;
  }
  strncpy((Char_t*)(lnam+6),"FSTR",4); lnum[6]=dig[2];
  strncpy((Char_t*)(lnam+7),"FSEN",4); lnum[7]=0;

  //divz
  if (dig[3]<1 || dig[3]>2) {
     ::Info("DtoM","Wrong z-division number : %d !",dig[3]);
     return 3;
  }
  strncpy((Char_t*)(lnam+8),"FSEZ",4); lnum[8]=dig[3];

  //divx
  if (dig[4]<1 || dig[4]>48) {
     ::Info("DtoM","Wrong x-division number : %d !",dig[4]);
     return 3;
  }
  strncpy((Char_t*)(lnam+9),"FSEX",4); lnum[9]=dig[4];
  strncpy((Char_t*)(lnam+10),"FPAD",4); lnum[10]=0;

  Gcvolu_t *gcvolu=geant3->Gcvolu(); gcvolu->nlevel=0;
  Int_t err=geant3->Glvolu(11,lnam,lnum);  //11-th level
  if (err) {
    ::Info("DtoM","Wrong volume !"); return 3;
  }
  Float_t l[3]={0.,0.,0}; geant3->Gdtom(l,g,1);
  return 0;
}


//_________________________________________________________________________
AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) {
  //
  //  The main constructor
  //
  if (InitGeo()) throw "AliTOFpidESD: can not initialize the geometry !\n";


  if (fTOFGeometry) delete fTOFGeometry;
  fTOFGeometry = new AliTOFGeometry();

  fR=378.; 
  fDy=fTOFGeometry->XPad(); fDz=fTOFGeometry->ZPad(); 
  fN=0; fEventN=0;

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

}


//_________________________________________________________________________

//_________________________________________________________________________
Int_t AliTOFpidESD::LoadClusters(const TFile *df) {

  //--------------------------------------------------------------------
  //This function loads the TOF clusters
  //--------------------------------------------------------------------

  if (!((TFile *)df)->IsOpen()) {
    Error("LoadClusters","file with the TOF digits has not been open !\n");
    return 1;
  }

  Char_t   name[100]; 
  sprintf(name,"TreeD%d",GetEventNumber());
  TTree *dTree=(TTree*)((TFile *)df)->Get(name);
  if (!dTree) { 
    Error("LoadClusters"," can't get the tree with digits !\n");
    return 1;
  }
  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();

    // fTOFGeometry->GetPos(dig,g);   // uncomment this 
    DtoM(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);
  }  

  delete dTree;
  return 0;
}

//_________________________________________________________________________
Int_t AliTOFpidESD::LoadClusters(TTree *dTree) {
  //--------------------------------------------------------------------
  //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();

    //fTOFGeometry->GetPos(dig,g);
    DtoM(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;

    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 dph=TMath::Abs(c->GetPhi()-phi);
      if (dph>TMath::Pi()) dph-=TMath::Pi();
      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;

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

    //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);
        continue;
      }
      p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(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"
	30
0f4a7374	31	#include <TGeant3.h> // For DtoM and InitGeo, should
	32	#include "../STRUCT/AliBODY.h" // be removed after we switch to
	33	#include "../STRUCT/AliFRAMEv2.h" // AliTOGeometry::GetPos
c630aafd	34	#include "AliTOFv2FHoles.h"
c630aafd	35
3f83f224	36	#include <stdlib.h>
3f83f224	37
c630aafd	38
	39	ClassImp(AliTOFpidESD)
	40
0f4a7374	41	//_________________________________________________________________________
c630aafd	42	static Int_t InitGeo() {
	43	//gSystem->Load("libgeant321");
	44	//new TGeant3("C++ Interface to Geant3");
	45
3f83f224	46	AliBODY *body = new AliBODY("BODY", "Alice envelop");
3f83f224	47	body->CreateGeometry();
c630aafd	48
3f83f224	49	AliFRAMEv2 *frame = new AliFRAMEv2("FRAME", "Space Frame");
	50	frame->SetHoles(1);
	51	frame->CreateGeometry();
c630aafd	52
3f83f224	53	AliTOF *tof = new AliTOFv2FHoles("TOF", "TOF with Holes");
3f83f224	54	tof->CreateGeometry();
c630aafd	55
	56	return 0;
	57	}
	58
	59	//_________________________________________________________________________
c630aafd	60	static Int_t DtoM(Int_t dig, Float_t g) {
3f83f224	61	const Int_t kMAX=13;
3f83f224	62	Int_t lnam[kMAX],lnum[kMAX];
c630aafd	63
c630aafd	64	TGeant3 geant3=(TGeant3)gMC;
3f83f224	65	if (!geant3) {::Info("DtoM","no geant3 found !"); return 1;}
c630aafd	66
	67	strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1;
	68	strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1;
	69
3f83f224	70	//4 padx if z<=0 then ...
0f4a7374	71	if ((dig[1]==4)\|\|(dig[1]==3)) dig[4]=AliTOFGeometry::NpadX()-dig[4];
3f83f224	72	else if (dig[1]==2) {
0f4a7374	73	if (dig[2]>7) dig[4]=AliTOFGeometry::NpadX()-dig[4];
3f83f224	74	else if (dig[2]==7) {
0f4a7374	75	if (dig[3]==1) dig[4]=AliTOFGeometry::NpadX()-dig[4];
3f83f224	76	else dig[4]+=1;
	77	} else dig[4]+=1;
	78	} else dig[4]+=1;
	79
	80	//3 padz
0f4a7374	81	if ((dig[1]==3)\|\|(dig[1]==4)) dig[3]=AliTOFGeometry::NpadZ()-dig[3];
3f83f224	82	else dig[3]+=1;
	83
	84	//2 strip
0f4a7374	85	if (dig[1]==0) dig[2]=AliTOFGeometry::NStripC()-dig[2];
0f4a7374	86	else if (dig[1]==1) dig[2]=AliTOFGeometry::NStripB()-dig[2];
3f83f224	87	else dig[2]+=1;
	88
	89	//1 module
0f4a7374	90	dig[1]=AliTOFGeometry::NPlates()-dig[1];
3f83f224	91
3f83f224	92	//0 sector
0f4a7374	93	dig[0]+=(AliTOFGeometry::NSectors()/2); dig[0]%=AliTOFGeometry::NSectors();
3f83f224	94	dig[0]+=1;
3f83f224	95
c630aafd	96	//sector
	97	switch (dig[0]) {
	98	case 1:
	99	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=9;
	100	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
	101	break;
	102	case 2:
	103	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=10;
	104	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
	105	break;
	106	case 3:
	107	strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=4;
	108	strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
	109	if (dig[1]>3) lnum[3]=2;
	110	break;
	111	case 4:
	112	strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=5;
	113	strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
	114	if (dig[1]>3) lnum[3]=2;
	115	break;
	116	case 5:
	117	strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=1;
	118	strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
	119	if (dig[1]>3) lnum[3]=2;
	120	break;
	121	case 6:
	122	strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=2;
	123	strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
	124	if (dig[1]>3) lnum[3]=2;
	125	break;
	126	case 7:
	127	strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=3;
	128	strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
	129	if (dig[1]>3) lnum[3]=2;
	130	break;
	131	case 8:
	132	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=1;
	133	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
	134	break;
	135	case 9:
	136	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=2;
	137	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
	138	break;
	139	case 10:
	140	strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=1;
	141	strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
	142	if (dig[1]>4) lnum[3]=2;
	143	break;
	144	case 11:
	145	strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=2;
	146	strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
	147	if (dig[1]>4) lnum[3]=2;
	148	break;
	149	case 12:
	150	strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=3;
	151	strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
	152	if (dig[1]>4) lnum[3]=2;
	153	break;
	154	case 13:
	155	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=3;
	156	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
	157	break;
	158	case 14:
	159	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=4;
160	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
161	break;
162	case 15:
163	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=5;
164	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
165	break;
166	case 16:
167	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=6;
168	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
169	break;
170	case 17:
171	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=7;
172	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
173	break;
174	case 18:
175	strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=8;
176	strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
177	break;
178	default:
3f83f224	179	::Info("DtoM","Wrong sector number : %d !",dig[0]);
c630aafd	180	return 2;
	181	}
	182	//module
	183	switch (dig[1]) {
	184	case 1:
	185	strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
	186	if (dig[0]==1 \|\| dig[0]==2 ) lnum[4]=2;
	187	if (dig[0]>=8 && dig[0]<=9) lnum[4]=2;
	188	if (dig[0]>=13 && dig[0]<=18) lnum[4]=2;
	189	strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
	190	break;
	191	case 2:
	192	strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
	193	if (dig[0]==1 \|\| dig[0]==2 ) lnum[4]=2;
	194	if (dig[0]>=8 && dig[0]<=9) lnum[4]=2;
	195	if (dig[0]>=13 && dig[0]<=18) lnum[4]=2;
	196	strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
	197	break;
	198	case 3:
	199	strncpy((Char_t*)(lnam+4),"FTOA",4); lnum[4]=0;
	200	strncpy((Char_t*)(lnam+5),"FLTA",4); lnum[5]=0;
	201	break;
	202	case 4:
	203	strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
	204	strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
	205	break;
	206	case 5:
	207	strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
	208	strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
	209	break;
	210	default:
3f83f224	211	::Info("DtoM","Wrong module number : %d !",dig[1]);
c630aafd	212	return 3;
	213	}
	214
	215	//strip
3f83f224	216	if (dig[2]<0/1/ \|\| dig[2]>19) {
3f83f224	217	::Info("DtoM","Wrong strip number : %d !",dig[2]);
c630aafd	218	return 3;
	219	}
	220	strncpy((Char_t*)(lnam+6),"FSTR",4); lnum[6]=dig[2];
	221	strncpy((Char_t*)(lnam+7),"FSEN",4); lnum[7]=0;
	222
	223	//divz
	224	if (dig[3]<1 \|\| dig[3]>2) {
3f83f224	225	::Info("DtoM","Wrong z-division number : %d !",dig[3]);
c630aafd	226	return 3;
	227	}
	228	strncpy((Char_t*)(lnam+8),"FSEZ",4); lnum[8]=dig[3];
	229
	230	//divx
	231	if (dig[4]<1 \|\| dig[4]>48) {
3f83f224	232	::Info("DtoM","Wrong x-division number : %d !",dig[4]);
c630aafd	233	return 3;
	234	}
	235	strncpy((Char_t*)(lnam+9),"FSEX",4); lnum[9]=dig[4];
	236	strncpy((Char_t*)(lnam+10),"FPAD",4); lnum[10]=0;
	237
	238	Gcvolu_t *gcvolu=geant3->Gcvolu(); gcvolu->nlevel=0;
	239	Int_t err=geant3->Glvolu(11,lnam,lnum); //11-th level
3f83f224	240	if (err) {
	241	::Info("DtoM","Wrong volume !"); return 3;
	242	}
c630aafd	243	Float_t l[3]={0.,0.,0}; geant3->Gdtom(l,g,1);
	244	return 0;
	245	}
	246
0f4a7374	247
	248	//_________________________________________________________________________
	249	AliTOFpidESD::AliTOFpidESD(Double_t param) throw (const Char_t ) {
	250	//
	251	// The main constructor
	252	//
	253	if (InitGeo()) throw "AliTOFpidESD: can not initialize the geometry !\n";
	254
	255
	256	if (fTOFGeometry) delete fTOFGeometry;
	257	fTOFGeometry = new AliTOFGeometry();
	258
	259	fR=378.;
	260	fDy=fTOFGeometry->XPad(); fDz=fTOFGeometry->ZPad();
	261	fN=0; fEventN=0;
	262
	263	fSigma=param[0];
	264	fRange=param[1];
	265
	266	}
	267
	268
	269	//_________________________________________________________________________
	270
	271	//_________________________________________________________________________
c630aafd	272	Int_t AliTOFpidESD::LoadClusters(const TFile *df) {
0f4a7374	273
c630aafd	274	//--------------------------------------------------------------------
	275	//This function loads the TOF clusters
	276	//--------------------------------------------------------------------
0f4a7374	277
c630aafd	278	if (!((TFile *)df)->IsOpen()) {
	279	Error("LoadClusters","file with the TOF digits has not been open !\n");
	280	return 1;
	281	}
	282
	283	Char_t name[100];
	284	sprintf(name,"TreeD%d",GetEventNumber());
	285	TTree dTree=(TTree)((TFile *)df)->Get(name);
	286	if (!dTree) {
	287	Error("LoadClusters"," can't get the tree with digits !\n");
	288	return 1;
	289	}
	290	TBranch *branch=dTree->GetBranch("TOF");
	291	if (!branch) {
	292	Error("LoadClusters"," can't get the branch with the TOF digits !\n");
	293	return 1;
	294	}
	295
	296	TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
	297	branch->SetAddress(&digits);
	298
	299	dTree->GetEvent(0);
	300	Int_t nd=digits->GetEntriesFast();
	301	Info("LoadClusters","number of digits: %d",nd);
	302
	303	for (Int_t i=0; i<nd; i++) {
	304	AliTOFdigit d=(AliTOFdigit)digits->UncheckedAt(i);
	305	Int_t dig[5]; Float_t g[3];
	306	dig[0]=d->GetSector();
	307	dig[1]=d->GetPlate();
	308	dig[2]=d->GetStrip();
	309	dig[3]=d->GetPadz();
	310	dig[4]=d->GetPadx();
	311
0f4a7374	312	// fTOFGeometry->GetPos(dig,g); // uncomment this
c630aafd	313	DtoM(dig,g);
	314
	315	Double_t h[5];
	316	h[0]=TMath::Sqrt(g[0]g[0]+g[1]g[1]);
	317	h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
	318	h[3]=d->GetTdc(); h[4]=d->GetAdc();
	319
3f83f224	320	AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
c630aafd	321	InsertCluster(cl);
	322	}
	323
	324	delete dTree;
	325	return 0;
	326	}
	327
0f4a7374	328	//_________________________________________________________________________
c630aafd	329	Int_t AliTOFpidESD::LoadClusters(TTree *dTree) {
	330	//--------------------------------------------------------------------
	331	//This function loads the TOF clusters
	332	//--------------------------------------------------------------------
	333	TBranch *branch=dTree->GetBranch("TOF");
	334	if (!branch) {
	335	Error("LoadClusters"," can't get the branch with the TOF digits !\n");
	336	return 1;
	337	}
	338
	339	TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
	340	branch->SetAddress(&digits);
	341
	342	dTree->GetEvent(0);
	343	Int_t nd=digits->GetEntriesFast();
	344	Info("LoadClusters","number of digits: %d",nd);
	345
	346	for (Int_t i=0; i<nd; i++) {
	347	AliTOFdigit d=(AliTOFdigit)digits->UncheckedAt(i);
	348	Int_t dig[5]; Float_t g[3];
	349	dig[0]=d->GetSector();
	350	dig[1]=d->GetPlate();
	351	dig[2]=d->GetStrip();
	352	dig[3]=d->GetPadz();
	353	dig[4]=d->GetPadx();
	354
0f4a7374	355	//fTOFGeometry->GetPos(dig,g);
c630aafd	356	DtoM(dig,g);
	357
	358	Double_t h[5];
	359	h[0]=TMath::Sqrt(g[0]g[0]+g[1]g[1]);
	360	h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
	361	h[3]=d->GetTdc(); h[4]=d->GetAdc();
	362
3f83f224	363	AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
c630aafd	364	InsertCluster(cl);
	365	}
	366
	367	return 0;
	368	}
	369
0f4a7374	370	//_________________________________________________________________________
c630aafd	371	void AliTOFpidESD::UnloadClusters() {
	372	//--------------------------------------------------------------------
	373	//This function unloads TOF clusters
	374	//--------------------------------------------------------------------
	375	for (Int_t i=0; i<fN; i++) delete fClusters[i];
	376	fN=0;
	377	}
	378
0f4a7374	379	//_________________________________________________________________________
c630aafd	380	Int_t AliTOFpidESD::InsertCluster(AliTOFcluster *c) {
	381	//--------------------------------------------------------------------
	382	//This function adds a cluster to the array of clusters sorted in Z
	383	//--------------------------------------------------------------------
	384	if (fN==kMaxCluster) {
	385	Error("InsertCluster","Too many clusters !\n");
	386	return 1;
	387	}
	388
	389	if (fN==0) {fClusters[fN++]=c; return 0;}
	390	Int_t i=FindClusterIndex(c->GetZ());
	391	memmove(fClusters+i+1 ,fClusters+i,(fN-i)sizeof(AliTOFcluster));
	392	fClusters[i]=c; fN++;
	393
	394	return 0;
	395	}
	396
0f4a7374	397	//_________________________________________________________________________
c630aafd	398	Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
	399	//--------------------------------------------------------------------
	400	// This function returns the index of the nearest cluster
	401	//--------------------------------------------------------------------
	402	if (fN==0) return 0;
	403	if (z <= fClusters[0]->GetZ()) return 0;
	404	if (z > fClusters[fN-1]->GetZ()) return fN;
	405	Int_t b=0, e=fN-1, m=(b+e)/2;
	406	for (; b<e; m=(b+e)/2) {
	407	if (z > fClusters[m]->GetZ()) b=m+1;
	408	else e=m;
	409	}
	410	return m;
	411	}
	412
3f83f224	413	static int cmp(const void p1, const void p2) {
	414	AliESDtrack t1=((AliESDtrack**)p1);
	415	AliESDtrack t2=((AliESDtrack**)p2);
	416	Double_t c1[15]; t1->GetExternalCovariance(c1);
	417	Double_t c2[15]; t2->GetExternalCovariance(c2);
	418	if (c1[0]c1[2] <c2[0]c2[2]) return -1;
	419	if (c1[0]c1[2]==c2[0]c2[2]) return 0;
	420	return 1;
	421	}
	422
c630aafd	423	//_________________________________________________________________________
	424	Int_t AliTOFpidESD::MakePID(AliESD *event)
	425	{
	426	//
	427	// This function calculates the "detector response" PID probabilities
	428	// Just for a bare hint...
	429
3f83f224	430	static const Double_t kMasses[]={
c630aafd	431	0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
	432	};
	433
	434	Int_t ntrk=event->GetNumberOfTracks();
3f83f224	435	AliESDtrack *tracks=new AliESDtrack[ntrk];
	436
	437	Int_t i;
	438	for (i=0; i<ntrk; i++) {
c630aafd	439	AliESDtrack *t=event->GetTrack(i);
3f83f224	440	tracks[i]=t;
	441	}
	442	qsort(tracks,ntrk,sizeof(AliESDtrack*),cmp);
	443
	444	Int_t nmatch=0;
	445	for (i=0; i<ntrk; i++) {
	446	AliESDtrack *t=tracks[i];
c630aafd	447
c630aafd	448	if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
c630aafd	449
	450	Double_t x,par[5]; t->GetExternalParameters(x,par);
	451	Double_t cov[15]; t->GetExternalCovariance(cov);
	452
3f83f224	453	Double_t dphi=(5TMath::Sqrt(cov[0]) + 0.5fDy + 2.5*TMath::Abs(par[2]))/fR;
3f83f224	454	Double_t dz=5TMath::Sqrt(cov[2]) + 0.5fDz + 2.5*TMath::Abs(par[3]);
c630aafd	455
	456	Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha();
	457	if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
	458	if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
	459	Double_t z=par[1];
	460
	461	Double_t d2max=1000.;
	462	Int_t index=-1;
	463	for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
	464	AliTOFcluster *c=fClusters[k];
	465
	466	if (c->GetZ() > z+dz) break;
	467	if (c->IsUsed()) continue;
	468
3f83f224	469	Double_t dph=TMath::Abs(c->GetPhi()-phi);
	470	if (dph>TMath::Pi()) dph-=TMath::Pi();
	471	if (dph>dphi) continue;
c630aafd	472
3f83f224	473	Double_t d2=dphdphfRfR + (c->GetZ()-z)(c->GetZ()-z);
c630aafd	474	if (d2 > d2max) continue;
	475
	476	d2max=d2;
	477	index=k;
	478	}
	479
	480	if (index<0) {
	481	//Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel()));
	482	continue;
	483	}
	484
	485	nmatch++;
	486
	487	AliTOFcluster *c=fClusters[index];
3f83f224	488	c->Use();
	489
	490	Double_t tof=50*c->GetTDC()+32; // in ps
	491	t->SetTOFsignal(tof);
	492	t->SetTOFcluster(c->GetIndex());
	493
	494	if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
c630aafd	495
c630aafd	496	Double_t time[10]; t->GetIntegratedTimes(time);
3f83f224	497
	498	//track length correction
	499	Double_t rc=TMath::Sqrt(c->GetR()c->GetR() + c->GetZ()c->GetZ());
	500	Double_t rt=TMath::Sqrt(xx + par[0]par[0] + par[1]*par[1]);
	501	Double_t dlt=rc-rt;
c630aafd	502
	503	Double_t p[10];
	504	Double_t mom=t->GetP();
	505	for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
3f83f224	506	Double_t mass=kMasses[j];
	507	Double_t dpp=0.01; //mean relative pt resolution;
	508	if (mom>0.5) dpp=0.01*mom;
	509	Double_t sigma=dpptime[j]/(1.+ mommom/(mass*mass));
c630aafd	510	sigma=TMath::Sqrt(sigmasigma + fSigmafSigma);
3f83f224	511
	512	time[j]+=dlt/3e-2TMath::Sqrt(mommom+mass*mass)/mom;
	513
c630aafd	514	if (TMath::Abs(tof-time[j]) > fRange*sigma) {
	515	p[j]=TMath::Exp(-0.5fRangefRange);
	516	continue;
	517	}
	518	p[j]=TMath::Exp(-0.5(tof-time[j])(tof-time[j])/(sigma*sigma));
	519	}
3f83f224	520	/*
	521	if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
	522	cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
	523	continue;
	524	}
	525	*/
c630aafd	526	t->SetTOFpid(p);
c630aafd	527
c630aafd	528	}
	529
	530	Info("MakePID","Number of matched ESD track: %d",nmatch);
	531
3f83f224	532	delete[] tracks;
3f83f224	533
c630aafd	534	return 0;
c630aafd	535	}
0f4a7374	536