[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliIDFFUtils.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.                  *
 **************************************************************************/

//Utils for identified fragmentation function (IDFF) analysis
//Author: Xianguo Lu (xianguo.lu@cern.ch)

#include "AliIDFFUtils.h"
#include "AliAODEvent.h"
#include "AliAODMCParticle.h"
#include "AliESDtrack.h"

#include "TClonesArray.h"

ClassImp(AliIDFFUtils);

AliPIDResponse * AliIDFFUtils::fPid=0x0;

Int_t AliIDFFUtils::PDG2Type(const Int_t pdg)
{
  //
  //conversion from pdg code to local definition of particle type
  //

  Int_t itype = kNOTSELECTED;
  switch(pdg){
  case 11:
    itype = kELECTRON;
    break;
  case 211:
    itype = kPION;
    break;
  case 2212:
    itype = kPROTON;
    break;
  case 321:
    itype = kKAON;
    break;
  default:
    itype = kNOTSELECTED;
    break;
  }
  return itype;
}

THnSparseD *AliIDFFUtils::GetTHn(const TString name)
{
  //
  //get THnSparseD
  //

  const Int_t nvar = 11;
  //                                   0       1              2              3               4     5      6       7       8         9     10   
  const TString atitle[nvar]={"TrackEta","JetPt", "TrackTPCsig", "Log10TrackP", "Log10TrackPt",  "z",  "xi",  "pdg",  "comb",   "tof",  "tpc"};
  const Int_t nbins[nvar]   ={         4,     15,          1200,          Nx(),             50,   30,    60,      7,      7,        7,      7};
  const Double_t xmins[nvar]={         0,      5,             0,        Xmin(),         Xmin(),    0,     0,   -3.5,   -3.5,     -3.5,   -3.5};
  const Double_t xmaxs[nvar]={       0.9,     20,           200,        Xmax(),         Xmax(),  1.2,     6,    3.5,    3.5,      3.5,    3.5};

  THnSparseD * hh = new THnSparseD(name,"", nvar, nbins, xmins, xmaxs);
  for(Int_t ia=0; ia<nvar; ia++){
    hh->GetAxis(ia)->SetTitle(atitle[ia]);
  }

  TAxis * ax = 0x0;
  Int_t nb = 0;

  //0: eta bin
  ax = hh->GetAxis(0);
  nb = ax->GetNbins();
  const Double_t etas[]={0, 0.2, 0.4, 0.6, 0.9};
  ax->Set(nb, etas);

  return hh;
}

/*
Bool_t AliIDFFUtils::HMPIDAcceptance(const AliAODTrack *track)
{
  //
  //check HMPID acceptance
  //From S. Pochybova
  //

  Double_t tEta = TMath::Abs(track->Eta());
  Double_t tPhi = track->Phi();
  if(tPhi < 0.) tPhi += TMath::TwoPi();
  if(tPhi > TMath::TwoPi()) tPhi -= TMath::TwoPi();

  if(tEta > 0.46){
    return kFALSE;
  }

  if(tPhi < 0.08 || tPhi > 1.12){
    return kFALSE;
  }

  return kTRUE;
}

Bool_t AliIDFFUtils::HMPIDQA(const AliAODTrack *track)
{
  //
  //check HMPID PID quality
  //From S. Pochybova
  //

  if(track->GetHMPIDsignal() <= 0.){
    return kFALSE;
  }
  
  //check track-quality cuts
  //dist_(mip-trk)
  //track variables
  Float_t tX, tY, tTh, tPh;
  //mip variables
  Float_t mpX, mpY;
  Int_t mpQ, mpNph;
  track->GetHMPIDtrk(tX, tY, tTh, tPh);
  track->GetHMPIDmip(mpX, mpY, mpQ, mpNph);
  const Double_t dist = TMath::Sqrt((tX-mpX)*(tX-mpX)+(tY-mpY)*(tY-mpY));

  //taking the pass 2 case for the moment
  if(dist > 1){
    //Printf("Track did not pass the distance cut");
    return kFALSE;
  }

  //cut on charge
  //have to check if this also varies with pass
  //taking the pass 2 case for the moment
  if(mpQ < 120){
    //Printf("Track did not pass the MipQ cut");
    return kFALSE;
  }

  return kTRUE;
}

Int_t AliIDFFUtils::HMPIDType(const AliAODTrack * track)
{
  //
  //return the (locally defined) particle type judged by HMPID
  //From S. Pochybova
  //

  if(!HMPIDAcceptance(track))
    return kNOTACCEPTED;

  if(!HMPIDQA(track))
    return kNOINFO;

  Double_t nsigma[]={-999,-999,-999};
  nsigma[kPION]     = fPid->NumberOfSigmasHMPID( track, AliPID::kPion);
  nsigma[kKAON]     = fPid->NumberOfSigmasHMPID( track, AliPID::kKaon);
  nsigma[kPROTON]   = fPid->NumberOfSigmasHMPID( track, AliPID::kProton);

  const Double_t inclusion=2;
  const Double_t exclusion=3;

  if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kPION;
  if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kKAON;
  if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion) return kPROTON;

  return kNOTSELECTED;
}
*/


Int_t AliIDFFUtils::TPCType(const AliAODTrack * trackptr)
{
  //
  //return the (locally defined) particle type judged by TPC
  //use tofmode for TPC mode
  //

  const AliPIDResponse::EDetPidStatus tpcstatus =  fPid->CheckPIDStatus(AliPIDResponse::kTPC, trackptr);
  if(tpcstatus != AliPIDResponse::kDetPidOk)
    return kNOINFO;

  Double_t nsigma[]={-999,-999,-999, -999};
  nsigma[kPION]     = fPid->NumberOfSigmasTPC( trackptr, AliPID::kPion);
  nsigma[kKAON]     = fPid->NumberOfSigmasTPC( trackptr, AliPID::kKaon);
  nsigma[kPROTON]   = fPid->NumberOfSigmasTPC( trackptr, AliPID::kProton);
  nsigma[kELECTRON] = fPid->NumberOfSigmasTPC( trackptr, AliPID::kElectron);

  //so that the effective region is really low momentum
  const Double_t inclusion=5;
  const Double_t exclusion=5;

  //don't destroy TPC signal shape below 120
  const Double_t maxsig = 150;
  if(trackptr->GetTPCsignal()> maxsig){
    if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kPION;
    if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kKAON;
    if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kPROTON;
    if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])<inclusion) return kELECTRON;
  }

  return kNOTSELECTED;
}

Int_t AliIDFFUtils::TOFType(const AliAODTrack * trackptr, const Int_t tofmode)
{
  //
  //return the (locally defined) particle type judged by TOF
  //

  //check kTOFout, kTIME, mismatch 
  const AliPIDResponse::EDetPidStatus tofstatus =  fPid->CheckPIDStatus(AliPIDResponse::kTOF, trackptr);
  if(tofstatus != AliPIDResponse::kDetPidOk)
    return kNOINFO;

  Double_t nsigma[]={-999,-999,-999, -999};
  nsigma[kPION]     = fPid->NumberOfSigmasTOF( trackptr, AliPID::kPion);
  nsigma[kKAON]     = fPid->NumberOfSigmasTOF( trackptr, AliPID::kKaon);
  nsigma[kPROTON]   = fPid->NumberOfSigmasTOF( trackptr, AliPID::kProton);
  nsigma[kELECTRON] = fPid->NumberOfSigmasTOF( trackptr, AliPID::kElectron);

  Double_t inclusion=-999;
  Double_t exclusion=-999;
  if(tofmode == 1){
    inclusion = 2;
    exclusion = 2;
  }
  else if(tofmode == 2){
    inclusion = 2;
    exclusion = 3;
  }
  else if(tofmode == 3){
    inclusion = 3;
    exclusion = 3;
  }
  else if(tofmode == 4){
    inclusion = 3;
    exclusion = 4;
  }
  else{
    printf("AliIDFFUtils::TOFType bad tofmode ! %d\n", tofmode); exit(1);
  }

  const Bool_t kpassEle = kTRUE;
  /*
  const Double_t cutEle = 1.5; 
  //tofmode = 1x then require electron exclusion cut
  if( tofmode == 4 ){
    if(TMath::Abs(nsigma[kELECTRON])> cutEle ){
      kpassEle = kTRUE;
    }
    else{
      kpassEle = kFALSE;
    }
  }
  */

  //cut on electron for pion because the precision of pion is good and the contamination of electron can not be ignored
  //+1 exclusion sigma in electron/pion to enforce better purity, otherwise not only pion, but also kaon is bias for jet pt 5-10 GeV/c
  if(TMath::Abs(nsigma[kPION])<inclusion     && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && kpassEle) return kPION;
  if(TMath::Abs(nsigma[kPION])>exclusion     && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion && kpassEle) return kKAON;
  if(TMath::Abs(nsigma[kPION])>exclusion     && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion && kpassEle) return kPROTON;
  if(TMath::Abs(nsigma[kPION])>exclusion+0.5 && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])<inclusion) return kELECTRON;

  return kNOTSELECTED;
}

Int_t AliIDFFUtils::CombineTPCTOF(const Int_t ktpc, const Int_t ktof)
{
  //tpc and tof, if <0 only noinfo or notselected
  if(ktpc == ktof)
    return ktpc;
  else if(ktpc < 0 && ktof >= 0 )
    return ktof;
  else if(ktof < 0 && ktpc >= 0)
    return ktpc;
  else
    return kNOTACCEPTED;
}

void AliIDFFUtils::FillTHn(THnSparseD * hh, Double_t jetpt, const AliAODTrack * track,  AliAODEvent *aodevent, const Int_t tofmode) //AliMCEvent * mcevent)
{
  //
  //fill variables
  //

  Int_t mcpdg = kNOINFO;

  TClonesArray *tca = dynamic_cast<TClonesArray*>(aodevent->FindListObject(AliAODMCParticle::StdBranchName()));
  if(tca){
    const Int_t mclabel = TMath::Abs(track->GetLabel());

    AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(mclabel));
    if(gentrack){
      //printf("null gentrack in AliIDFFUtils::FillTHn mclabel %d\n", mclabel);
      mcpdg = PDG2Type(TMath::Abs(gentrack->GetPdgCode()));
    }
  }
  
  //===========================================================================================================

  //use tofmode for tpcmode
  const Int_t ktpc = TPCType(track);
  const Int_t ktof = TOFType(track, tofmode);

  //fake kcomb by pretpc+pretof
  const Int_t kcomb = CombineTPCTOF(ktpc, ktof);

  //const Int_t khmpid = HMPIDType(track);

  //===========================================================================================================

  const Double_t eps = 1e-6;
  const Double_t tracketa = TMath::Abs(track->Eta());
  Double_t tracktpc = track->GetTPCsignal();
  if(tracktpc>=200)
    tracktpc = 200 - eps;
  const Double_t tracklogmom = TMath::Log10(track->P());
  const Double_t tracklogpt  = TMath::Log10(track->Pt());

  Double_t zz = -999;
  Double_t xi = -999;
  if(jetpt<-990){
    jetpt = zz = xi = eps;
  }
  //from Oliver
  else if(track->Pt()>(1-eps)*jetpt && track->Pt()<(1+eps)*jetpt){ // case z=1 : move entry to last histo bin <1
    zz  = 1-eps;
    xi = eps;
  }
  else if(jetpt>0){
    zz = track->Pt()/jetpt;
    xi = TMath::Log(1/zz);
  }

  //===========================================================================================================

  const Double_t vars[]={tracketa, jetpt, tracktpc, tracklogmom, tracklogpt, zz, xi, (Double_t)mcpdg, (Double_t)kcomb, (Double_t)ktof, (Double_t)ktpc};
  
  hh->Fill(vars);
}

//_________________________________________________________________________________
Bool_t AliIDFFUtils::TPCCutPIDN(const AliAODTrack * track)
{
  //
  //TPC Cut kPIDN
  //

  //assuming this cut is particle type independent, then it is fine
  //need further investigation
  if(track->GetTPCsignalN()<60){
    return kFALSE;
  }

  return kTRUE;
}

//_________________________________________________________________________________
Bool_t AliIDFFUtils::TPCCutMIGeo(const AliAODTrack * track, const AliVEvent* evt, TTreeStream * streamer)
{
  //
  //TPC Cut MIGeo
  //

  Short_t sign = track->Charge();
  Double_t xyz[50];
  Double_t pxpypz[50];
  Double_t cv[100];

  track->GetXYZ(xyz);
  track->GetPxPyPz(pxpypz);

  AliExternalTrackParam * par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
  AliESDtrack dummy;

  Double_t varGeom = dummy.GetLengthInActiveZone(par,3,236, evt->GetMagneticField(), 0,0);
  Double_t varNcr  = track->GetTPCClusterInfo(3,1);
  Double_t varNcls = track->GetTPCsignalN();

  Bool_t cutGeom     = varGeom > 1.*(130-5*TMath::Abs(1./track->Pt()));
  Bool_t cutNcr      = varNcr  > 0.85*(130-5*TMath::Abs(1./track->Pt()));
  Bool_t cutNcls     = varNcls > 0.7*(130-5*TMath::Abs(1./track->Pt()));

  Bool_t kout = cutGeom && cutNcr && cutNcls;

  if(streamer){
    Double_t dedx = track->GetTPCsignal();

    (*streamer)<<"tree"<<
      "param.="<< par<<

      "varGeom="<<varGeom<<
      "varNcr="<<varNcr<<
      "varNcls="<<varNcls<<

      "cutGeom="<<cutGeom<<
      "cutNcr="<<cutNcr<<
      "cutNcls="<<cutNcls<<

      "kout="<<kout<<
      "dedx="<<dedx<<

      "\n";
  }

  delete par;

  return kout;
}
Commit	Line	Data
7f0c28ff ML	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	//Utils for identified fragmentation function (IDFF) analysis
	17	//Author: Xianguo Lu (xianguo.lu@cern.ch)
	18
	19	#include "AliIDFFUtils.h"
	20	#include "AliAODEvent.h"
	21	#include "AliAODMCParticle.h"
	22	#include "AliESDtrack.h"
	23
	24	#include "TClonesArray.h"
	25
	26	ClassImp(AliIDFFUtils);
	27
	28	AliPIDResponse * AliIDFFUtils::fPid=0x0;
	29
	30	Int_t AliIDFFUtils::PDG2Type(const Int_t pdg)
	31	{
	32	//
	33	//conversion from pdg code to local definition of particle type
	34	//
	35
	36	Int_t itype = kNOTSELECTED;
	37	switch(pdg){
	38	case 11:
	39	itype = kELECTRON;
	40	break;
	41	case 211:
	42	itype = kPION;
	43	break;
	44	case 2212:
	45	itype = kPROTON;
	46	break;
	47	case 321:
	48	itype = kKAON;
	49	break;
	50	default:
	51	itype = kNOTSELECTED;
	52	break;
	53	}
	54	return itype;
	55	}
	56
	57	THnSparseD *AliIDFFUtils::GetTHn(const TString name)
	58	{
	59	//
	60	//get THnSparseD
	61	//
	62
1bb86680	63	const Int_t nvar = 11;
	64	// 0 1 2 3 4 5 6 7 8 9 10
	65	const TString atitle[nvar]={"TrackEta","JetPt", "TrackTPCsig", "Log10TrackP", "Log10TrackPt", "z", "xi", "pdg", "comb", "tof", "tpc"};
	66	const Int_t nbins[nvar] ={ 4, 15, 1200, Nx(), 50, 30, 60, 7, 7, 7, 7};
	67	const Double_t xmins[nvar]={ 0, 5, 0, Xmin(), Xmin(), 0, 0, -3.5, -3.5, -3.5, -3.5};
	68	const Double_t xmaxs[nvar]={ 0.9, 20, 200, Xmax(), Xmax(), 1.2, 6, 3.5, 3.5, 3.5, 3.5};
7f0c28ff ML	69
	70	THnSparseD * hh = new THnSparseD(name,"", nvar, nbins, xmins, xmaxs);
	71	for(Int_t ia=0; ia<nvar; ia++){
	72	hh->GetAxis(ia)->SetTitle(atitle[ia]);
	73	}
	74
	75	TAxis * ax = 0x0;
	76	Int_t nb = 0;
	77
	78	//0: eta bin
	79	ax = hh->GetAxis(0);
	80	nb = ax->GetNbins();
	81	const Double_t etas[]={0, 0.2, 0.4, 0.6, 0.9};
	82	ax->Set(nb, etas);
	83
	84	return hh;
	85	}
	86
1bb86680	87	/*
7f0c28ff ML	88	Bool_t AliIDFFUtils::HMPIDAcceptance(const AliAODTrack *track)
	89	{
	90	//
	91	//check HMPID acceptance
	92	//From S. Pochybova
	93	//
	94
	95	Double_t tEta = TMath::Abs(track->Eta());
	96	Double_t tPhi = track->Phi();
	97	if(tPhi < 0.) tPhi += TMath::TwoPi();
	98	if(tPhi > TMath::TwoPi()) tPhi -= TMath::TwoPi();
	99
	100	if(tEta > 0.46){
	101	return kFALSE;
	102	}
	103
	104	if(tPhi < 0.08 \|\| tPhi > 1.12){
	105	return kFALSE;
	106	}
	107
	108	return kTRUE;
	109	}
	110
	111	Bool_t AliIDFFUtils::HMPIDQA(const AliAODTrack *track)
	112	{
	113	//
	114	//check HMPID PID quality
	115	//From S. Pochybova
	116	//
	117
	118	if(track->GetHMPIDsignal() <= 0.){
	119	return kFALSE;
	120	}
	121
	122	//check track-quality cuts
	123	//dist_(mip-trk)
	124	//track variables
	125	Float_t tX, tY, tTh, tPh;
	126	//mip variables
	127	Float_t mpX, mpY;
	128	Int_t mpQ, mpNph;
	129	track->GetHMPIDtrk(tX, tY, tTh, tPh);
	130	track->GetHMPIDmip(mpX, mpY, mpQ, mpNph);
	131	const Double_t dist = TMath::Sqrt((tX-mpX)(tX-mpX)+(tY-mpY)(tY-mpY));
	132
	133	//taking the pass 2 case for the moment
	134	if(dist > 1){
	135	//Printf("Track did not pass the distance cut");
	136	return kFALSE;
	137	}
	138
	139	//cut on charge
	140	//have to check if this also varies with pass
	141	//taking the pass 2 case for the moment
	142	if(mpQ < 120){
	143	//Printf("Track did not pass the MipQ cut");
	144	return kFALSE;
	145	}
	146
	147	return kTRUE;
	148	}
	149
	150	Int_t AliIDFFUtils::HMPIDType(const AliAODTrack * track)
	151	{
152	//
153	//return the (locally defined) particle type judged by HMPID
154	//From S. Pochybova
155	//
156
157	if(!HMPIDAcceptance(track))
158	return kNOTACCEPTED;
159
160	if(!HMPIDQA(track))
161	return kNOINFO;
162
163	Double_t nsigma[]={-999,-999,-999};
164	nsigma[kPION] = fPid->NumberOfSigmasHMPID( track, AliPID::kPion);
165	nsigma[kKAON] = fPid->NumberOfSigmasHMPID( track, AliPID::kKaon);
166	nsigma[kPROTON] = fPid->NumberOfSigmasHMPID( track, AliPID::kProton);
167
168	const Double_t inclusion=2;
169	const Double_t exclusion=3;
170
171	if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kPION;
172	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kKAON;
173	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion) return kPROTON;
174
175	return kNOTSELECTED;
176	}
1bb86680	177	*/
	178
	179
	180	Int_t AliIDFFUtils::TPCType(const AliAODTrack * trackptr)
	181	{
	182	//
	183	//return the (locally defined) particle type judged by TPC
	184	//use tofmode for TPC mode
	185	//
	186
	187	const AliPIDResponse::EDetPidStatus tpcstatus = fPid->CheckPIDStatus(AliPIDResponse::kTPC, trackptr);
	188	if(tpcstatus != AliPIDResponse::kDetPidOk)
	189	return kNOINFO;
	190
	191	Double_t nsigma[]={-999,-999,-999, -999};
	192	nsigma[kPION] = fPid->NumberOfSigmasTPC( trackptr, AliPID::kPion);
	193	nsigma[kKAON] = fPid->NumberOfSigmasTPC( trackptr, AliPID::kKaon);
	194	nsigma[kPROTON] = fPid->NumberOfSigmasTPC( trackptr, AliPID::kProton);
	195	nsigma[kELECTRON] = fPid->NumberOfSigmasTPC( trackptr, AliPID::kElectron);
	196
	197	//so that the effective region is really low momentum
	198	const Double_t inclusion=5;
	199	const Double_t exclusion=5;
	200
	201	//don't destroy TPC signal shape below 120
	202	const Double_t maxsig = 150;
	203	if(trackptr->GetTPCsignal()> maxsig){
	204	if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kPION;
	205	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kKAON;
	206	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kPROTON;
	207	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])<inclusion) return kELECTRON;
	208	}
	209
	210	return kNOTSELECTED;
	211	}
7f0c28ff ML	212
	213	Int_t AliIDFFUtils::TOFType(const AliAODTrack * trackptr, const Int_t tofmode)
	214	{
	215	//
	216	//return the (locally defined) particle type judged by TOF
	217	//
	218
	219	//check kTOFout, kTIME, mismatch
	220	const AliPIDResponse::EDetPidStatus tofstatus = fPid->CheckPIDStatus(AliPIDResponse::kTOF, trackptr);
	221	if(tofstatus != AliPIDResponse::kDetPidOk)
	222	return kNOINFO;
	223
1bb86680	224	Double_t nsigma[]={-999,-999,-999, -999};
7f0c28ff ML	225	nsigma[kPION] = fPid->NumberOfSigmasTOF( trackptr, AliPID::kPion);
	226	nsigma[kKAON] = fPid->NumberOfSigmasTOF( trackptr, AliPID::kKaon);
	227	nsigma[kPROTON] = fPid->NumberOfSigmasTOF( trackptr, AliPID::kProton);
1bb86680	228	nsigma[kELECTRON] = fPid->NumberOfSigmasTOF( trackptr, AliPID::kElectron);
7f0c28ff ML	229
	230	Double_t inclusion=-999;
	231	Double_t exclusion=-999;
1bb86680	232	if(tofmode == 1){
	233	inclusion = 2;
	234	exclusion = 2;
7f0c28ff	235	}
1bb86680	236	else if(tofmode == 2){
7f0c28ff ML	237	inclusion = 2;
	238	exclusion = 3;
	239	}
1bb86680	240	else if(tofmode == 3){
1bb86680	241	inclusion = 3;
7f0c28ff ML	242	exclusion = 3;
7f0c28ff ML	243	}
1bb86680	244	else if(tofmode == 4){
	245	inclusion = 3;
	246	exclusion = 4;
	247	}
7f0c28ff ML	248	else{
	249	printf("AliIDFFUtils::TOFType bad tofmode ! %d\n", tofmode); exit(1);
	250	}
	251
1bb86680	252	const Bool_t kpassEle = kTRUE;
	253	/*
	254	const Double_t cutEle = 1.5;
	255	//tofmode = 1x then require electron exclusion cut
	256	if( tofmode == 4 ){
	257	if(TMath::Abs(nsigma[kELECTRON])> cutEle ){
	258	kpassEle = kTRUE;
	259	}
	260	else{
	261	kpassEle = kFALSE;
	262	}
	263	}
	264	*/
	265
	266	//cut on electron for pion because the precision of pion is good and the contamination of electron can not be ignored
	267	//+1 exclusion sigma in electron/pion to enforce better purity, otherwise not only pion, but also kaon is bias for jet pt 5-10 GeV/c
	268	if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && kpassEle) return kPION;
	269	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion && kpassEle) return kKAON;
	270	if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion && kpassEle) return kPROTON;
	271	if(TMath::Abs(nsigma[kPION])>exclusion+0.5 && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])<inclusion) return kELECTRON;
7f0c28ff ML	272
	273	return kNOTSELECTED;
	274	}
	275
1bb86680	276	Int_t AliIDFFUtils::CombineTPCTOF(const Int_t ktpc, const Int_t ktof)
	277	{
	278	//tpc and tof, if <0 only noinfo or notselected
	279	if(ktpc == ktof)
	280	return ktpc;
	281	else if(ktpc < 0 && ktof >= 0 )
	282	return ktof;
	283	else if(ktof < 0 && ktpc >= 0)
	284	return ktpc;
	285	else
	286	return kNOTACCEPTED;
	287	}
7f0c28ff ML	288
	289	void AliIDFFUtils::FillTHn(THnSparseD * hh, Double_t jetpt, const AliAODTrack * track, AliAODEvent aodevent, const Int_t tofmode) //AliMCEvent mcevent)
	290	{
	291	//
	292	//fill variables
	293	//
	294
	295	Int_t mcpdg = kNOINFO;
	296
	297	TClonesArray tca = dynamic_cast<TClonesArray>(aodevent->FindListObject(AliAODMCParticle::StdBranchName()));
	298	if(tca){
	299	const Int_t mclabel = TMath::Abs(track->GetLabel());
	300
	301	AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(mclabel));
	302	if(gentrack){
	303	//printf("null gentrack in AliIDFFUtils::FillTHn mclabel %d\n", mclabel);
	304	mcpdg = PDG2Type(TMath::Abs(gentrack->GetPdgCode()));
	305	}
	306	}
	307
	308	//===========================================================================================================
	309
1bb86680	310	//use tofmode for tpcmode
1bb86680	311	const Int_t ktpc = TPCType(track);
7f0c28ff	312	const Int_t ktof = TOFType(track, tofmode);
1bb86680	313
	314	//fake kcomb by pretpc+pretof
	315	const Int_t kcomb = CombineTPCTOF(ktpc, ktof);
	316
	317	//const Int_t khmpid = HMPIDType(track);
7f0c28ff ML	318
	319	//===========================================================================================================
	320
	321	const Double_t eps = 1e-6;
	322	const Double_t tracketa = TMath::Abs(track->Eta());
	323	Double_t tracktpc = track->GetTPCsignal();
	324	if(tracktpc>=200)
	325	tracktpc = 200 - eps;
	326	const Double_t tracklogmom = TMath::Log10(track->P());
	327	const Double_t tracklogpt = TMath::Log10(track->Pt());
	328
	329	Double_t zz = -999;
	330	Double_t xi = -999;
	331	if(jetpt<-990){
	332	jetpt = zz = xi = eps;
	333	}
	334	//from Oliver
	335	else if(track->Pt()>(1-eps)jetpt && track->Pt()<(1+eps)jetpt){ // case z=1 : move entry to last histo bin <1
	336	zz = 1-eps;
	337	xi = eps;
	338	}
	339	else if(jetpt>0){
	340	zz = track->Pt()/jetpt;
	341	xi = TMath::Log(1/zz);
	342	}
	343
	344	//===========================================================================================================
	345
1bb86680	346	const Double_t vars[]={tracketa, jetpt, tracktpc, tracklogmom, tracklogpt, zz, xi, (Double_t)mcpdg, (Double_t)kcomb, (Double_t)ktof, (Double_t)ktpc};
7f0c28ff ML	347
	348	hh->Fill(vars);
	349	}
	350
	351	//_________________________________________________________________________________
	352	Bool_t AliIDFFUtils::TPCCutPIDN(const AliAODTrack * track)
	353	{
	354	//
	355	//TPC Cut kPIDN
	356	//
	357
	358	//assuming this cut is particle type independent, then it is fine
	359	//need further investigation
	360	if(track->GetTPCsignalN()<60){
	361	return kFALSE;
	362	}
	363
	364	return kTRUE;
	365	}
	366
	367	//_________________________________________________________________________________
	368	Bool_t AliIDFFUtils::TPCCutMIGeo(const AliAODTrack * track, const AliVEvent* evt, TTreeStream * streamer)
	369	{
	370	//
	371	//TPC Cut MIGeo
	372	//
	373
	374	Short_t sign = track->Charge();
	375	Double_t xyz[50];
	376	Double_t pxpypz[50];
	377	Double_t cv[100];
	378
	379	track->GetXYZ(xyz);
	380	track->GetPxPyPz(pxpypz);
	381
	382	AliExternalTrackParam * par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
	383	AliESDtrack dummy;
	384
	385	Double_t varGeom = dummy.GetLengthInActiveZone(par,3,236, evt->GetMagneticField(), 0,0);
	386	Double_t varNcr = track->GetTPCClusterInfo(3,1);
	387	Double_t varNcls = track->GetTPCsignalN();
	388
	389	Bool_t cutGeom = varGeom > 1.(130-5TMath::Abs(1./track->Pt()));
	390	Bool_t cutNcr = varNcr > 0.85(130-5TMath::Abs(1./track->Pt()));
	391	Bool_t cutNcls = varNcls > 0.7(130-5TMath::Abs(1./track->Pt()));
	392
	393	Bool_t kout = cutGeom && cutNcr && cutNcls;
	394
	395	if(streamer){
	396	Double_t dedx = track->GetTPCsignal();
	397
	398	(*streamer)<<"tree"<<
	399	"param.="<< par<<
	400
	401	"varGeom="<<varGeom<<
	402	"varNcr="<<varNcr<<
	403	"varNcls="<<varNcls<<
	404
	405	"cutGeom="<<cutGeom<<
	406	"cutNcr="<<cutNcr<<
	407	"cutNcls="<<cutNcls<<
	408
	409	"kout="<<kout<<
	410	"dedx="<<dedx<<
411
412	"\n";
413	}
414
415	delete par;
416
417	return kout;
418	}
419