[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidEMCAL.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.                  *
**************************************************************************/
//
// Class for EMCAL PID
// Implements the abstract base class AliHFEpidBase
// IsInitialized() does the PID decision with energy and 
// momentum matching (e/p) 
// 
// Authors:
//  Shingo Sakai 
//   
//
//
#include <TMath.h>
#include "AliESDInputHandler.h"
#include "AliESDpid.h"

#include "AliHFEdetPIDqa.h"
#include "AliHFEpidEMCAL.h"
#include "AliHFEpidQAmanager.h"

#include "AliHFEemcalPIDqa.h"
//#include "AliVCluster.h"
//#include "AliVCaloCells.h"
//#include "AliVEvent.h"
#include "AliLog.h"
#include "AliPID.h"
//#include "AliESDEvent.h"
//#include "AliESDtrack.h"
#include "AliHFEpidEMCAL.h"

ClassImp(AliHFEpidEMCAL)

//___________________________________________________________________
AliHFEpidEMCAL::AliHFEpidEMCAL():
  AliHFEpidBase()
  , fPID(NULL)
  , feopMim(0.8)
  , feopMax(1.4)
{
  //
  // Constructor
  //
} 

//___________________________________________________________________
AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
  AliHFEpidBase(name)
  , fPID(NULL)
  , feopMim(0.8)
  , feopMax(1.4)
{
  //
  // Constructor
  //
}

//___________________________________________________________________
AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
  AliHFEpidBase("")
  , fPID(NULL)
  , feopMim(0.8)
  , feopMax(1.4)
{  
  // 
  // Copy operator
  //

  c.Copy(*this);
}
//___________________________________________________________________
AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
  //
  // Assignment operator
  //

  if(this != &ref){
    ref.Copy(*this);
  }

  return *this;
}
//___________________________________________________________________
AliHFEpidEMCAL::~AliHFEpidEMCAL(){
  //
  // Destructor
  //
  if(fPID) delete fPID;
}
//___________________________________________________________________
void AliHFEpidEMCAL::Copy(TObject &ref) const {
  //
  // Performs the copying of the object
  //
  AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);

  target.fPID = fPID;          
  target.feopMax = feopMax;
  target.feopMim = feopMim;

  AliHFEpidBase::Copy(ref);
}
//___________________________________________________________________
Bool_t AliHFEpidEMCAL::InitializePID(Int_t /*run*/){
  //
  // InitializePID: EMCAL experts have to implement code here
  //
  return kTRUE;
}

//___________________________________________________________________
Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
{ //Function to a return a code indicating whether or not an electron candidate is selected
  //
  //
  if(track==NULL)return 0;

  if(!fkPIDResponse) return 0;
  AliDebug(2, "PID object available");
  // EMCal not fESDpid  (s.s Feb. 11)
	
  const AliVTrack *trk = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
  if (trk == NULL) return 0;
  //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
  if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
  AliDebug(2, "Track Has EMCAL PID");
  
  //if(pidqa) 
  pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
  // not QA for EMCal, will be added (s.s Feb. 11)


     double feopMimCut = 0.0; 
     double feopMaxCut = 0.0; 

     feopMimCut = feopMim;
     feopMaxCut = feopMax;

     if(feopMax >900)  // not use fix cut
       {
        feopMimCut = CalEopCutMim(track->GetRecTrack(),0);
        feopMaxCut = CalEopCutMax(track->GetRecTrack(),0); 
       }
     if(feopMax < -900)  // not use fix cut for MC
       {
        feopMimCut = CalEopCutMim(track->GetRecTrack(),1);
        feopMaxCut = CalEopCutMax(track->GetRecTrack(),1); 
       }

  //printf("eop cuts org; %g; %g \n",feopMim,feopMax);
  //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);

  Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
  AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
  AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax));
  Int_t pdg = 0;
  //if(eop>feopMim && eop<feopMax){
  //if(eop>feopMim && eop<feopMax){
  if(eop>feopMimCut && eop<feopMaxCut){
    pdg = 11;
    //printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop);
    //if(pidqa) 
    pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
  }
  else
  {
   pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
  } 

    AliDebug(1, Form("eID %g ; %d \n",eop,pdg));  
    //printf("eID %g ; %d \n",eop,pdg);  

  return pdg;
}


//__________________________________________________________________________
Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const
{
  double maxCut = 0.0;
  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack==NULL)return maxCut;
  double Pt = esdtrack->Pt();

  if(flageop<0.5)
    {
     double meanP[3] = {0.991,1.0819,0.235}; 
     double sigP[3] = {6.52e-02,2.04e-01,3.34e-01}; 
     double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt); 
     double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt); 
     maxCut = mean+3.0*sig; 
    }
  else
    {
     double meanP[3] = {0.99,1.299,0.35}; 
     double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; 
     double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt); 
     double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt); 
     maxCut = mean+3.0*sig; 
    }
     return maxCut;
}

Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
{
  double mimCut = 0.0;
  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack==NULL)return mimCut;
  double Pt = esdtrack->Pt();

  if(flageop<0.5) // real
     { 
      //printf("real"); 
     double meanP[3] = {0.991,1.0819,0.235}; 
     double sigP[3] = {6.52e-02,2.04e-01,3.34e-01}; 
     double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt); 
     double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt); 
     mimCut = mean-2.0*sig;
    }
   else // MC
    { 
      //printf("MC"); 
     double meanP[3] = {0.99,1.299,0.35}; 
     double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; 
     double mean = meanP[0]*tanh(meanP[1]+meanP[2]*Pt); 
     double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt); 
     mimCut = mean-2.0*sig;
    }
     return mimCut;
}

//___________________________________________________________________________
Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
{ // Returns e/p & Rmatch

  Double_t matchclsE = 9999.9;
  double feop = -9999.9;

  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  if(esdtrack==NULL)return feop;
  const AliESDEvent *evt = esdtrack->GetESDEvent();

   Int_t icl = esdtrack->GetEMCALcluster();

   AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
   if(!cluster->IsEMCAL()) {return feop;}

       matchclsE = cluster->E();
  

   if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();

   return feop;

}


/*
//____________________________________________________________________________________
Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
{ // Returns e/p if an electron is matched

  Float_t  clsPos[3];
  Double_t trkPos[3];
  Double_t matchclsE = 9999.9;

  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
  AliESDEvent *evt = esdtrack->GetESDEvent();
  Double_t  magF = evt->GetMagneticField();
  Double_t magSign = 1.0;
  if(magF<0)magSign = -1.0;
  //printf("magF ; %g ; %g \n", magF,magSign);
 
  if (!TGeoGlobalMagField::Instance()->GetField()) {
	  printf("Loading field map...\n");
	  //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
	  AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
	  TGeoGlobalMagField::Instance()->SetField(field);
	            }

  AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
  Double_t fieldB[3]; 
  emctrack->GetBxByBz(fieldB);
  //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);

  for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
   {

   AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
   if(!cluster->IsEMCAL()) continue;

   cluster->GetPosition(clsPos);
   if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  continue;
   emctrack->GetXYZ(trkPos);

   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);

   Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
   Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching

   double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));

   if(rmatch<0.02)
      {
       matchclsE = cluster->E();
      }
  }
   delete emctrack;

   double feop = -9999.9;
   if(matchclsE<9999) feop = matchclsE/esdtrack->P();

   //   if(feop!=-9999.9)printf("%f\n",feop) ; 
   return feop;

}
*/
Commit	Line	Data
c2690925	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	// Class for EMCAL PID
	17	// Implements the abstract base class AliHFEpidBase
	18	// IsInitialized() does the PID decision with energy and
	19	// momentum matching (e/p)
	20	//
	21	// Authors:
	22	// Shingo Sakai
	23	//
	24	//
	25	//
	26	#include <TMath.h>
	27	#include "AliESDInputHandler.h"
c2690925	28	#include "AliESDpid.h"
	29
	30	#include "AliHFEdetPIDqa.h"
e156c3bb	31	#include "AliHFEpidEMCAL.h"
c2690925	32	#include "AliHFEpidQAmanager.h"
c2690925	33
e156c3bb	34	#include "AliHFEemcalPIDqa.h"
	35	//#include "AliVCluster.h"
	36	//#include "AliVCaloCells.h"
	37	//#include "AliVEvent.h"
c2690925	38	#include "AliLog.h"
c2690925	39	#include "AliPID.h"
e156c3bb	40	//#include "AliESDEvent.h"
e156c3bb	41	//#include "AliESDtrack.h"
c2690925	42	#include "AliHFEpidEMCAL.h"
	43
	44	ClassImp(AliHFEpidEMCAL)
	45
	46	//___________________________________________________________________
	47	AliHFEpidEMCAL::AliHFEpidEMCAL():
	48	AliHFEpidBase()
	49	, fPID(NULL)
	50	, feopMim(0.8)
	51	, feopMax(1.4)
	52	{
	53	//
	54	// Constructor
	55	//
	56	}
	57
	58	//___________________________________________________________________
	59	AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
	60	AliHFEpidBase(name)
	61	, fPID(NULL)
	62	, feopMim(0.8)
	63	, feopMax(1.4)
	64	{
	65	//
	66	// Constructor
	67	//
	68	}
	69
	70	//___________________________________________________________________
	71	AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
	72	AliHFEpidBase("")
	73	, fPID(NULL)
	74	, feopMim(0.8)
	75	, feopMax(1.4)
	76	{
	77	//
	78	// Copy operator
	79	//
	80
	81	c.Copy(*this);
	82	}
	83	//___________________________________________________________________
	84	AliHFEpidEMCAL &AliHFEpidEMCAL::operator=(const AliHFEpidEMCAL &ref){
	85	//
	86	// Assignment operator
	87	//
	88
	89	if(this != &ref){
	90	ref.Copy(*this);
	91	}
	92
	93	return *this;
	94	}
	95	//___________________________________________________________________
	96	AliHFEpidEMCAL::~AliHFEpidEMCAL(){
	97	//
	98	// Destructor
	99	//
	100	if(fPID) delete fPID;
	101	}
	102	//___________________________________________________________________
	103	void AliHFEpidEMCAL::Copy(TObject &ref) const {
	104	//
	105	// Performs the copying of the object
106	//
107	AliHFEpidEMCAL &target = dynamic_cast<AliHFEpidEMCAL &>(ref);
108
109	target.fPID = fPID;
110	target.feopMax = feopMax;
111	target.feopMim = feopMim;
112
113	AliHFEpidBase::Copy(ref);
114	}
115	//___________________________________________________________________
8c1c76e9	116	Bool_t AliHFEpidEMCAL::InitializePID(Int_t /run/){
c2690925	117	//
	118	// InitializePID: EMCAL experts have to implement code here
	119	//
	120	return kTRUE;
	121	}
	122
	123	//___________________________________________________________________
e156c3bb	124	Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject track, AliHFEpidQAmanager pidqa) const
c2690925	125	{ //Function to a return a code indicating whether or not an electron candidate is selected
	126	//
	127	//
63bffcf1	128	if(track==NULL)return 0;
63bffcf1	129
8c1c76e9	130	if(!fkPIDResponse) return 0;
c2690925	131	AliDebug(2, "PID object available");
c2690925	132	// EMCal not fESDpid (s.s Feb. 11)
63bffcf1	133
	134	const AliVTrack trk = dynamic_cast<const AliVTrack >(track->GetRecTrack());
	135	if (trk == NULL) return 0;
e156c3bb	136	//AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
63bffcf1	137	if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
c2690925	138	AliDebug(2, "Track Has EMCAL PID");
e156c3bb	139
	140	//if(pidqa)
	141	pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
c2690925	142	// not QA for EMCal, will be added (s.s Feb. 11)
c2690925	143
8c1c76e9	144
	145	double feopMimCut = 0.0;
	146	double feopMaxCut = 0.0;
	147
	148	feopMimCut = feopMim;
	149	feopMaxCut = feopMax;
	150
	151	if(feopMax >900) // not use fix cut
	152	{
	153	feopMimCut = CalEopCutMim(track->GetRecTrack(),0);
	154	feopMaxCut = CalEopCutMax(track->GetRecTrack(),0);
	155	}
	156	if(feopMax < -900) // not use fix cut for MC
	157	{
	158	feopMimCut = CalEopCutMim(track->GetRecTrack(),1);
	159	feopMaxCut = CalEopCutMax(track->GetRecTrack(),1);
	160	}
	161
	162	//printf("eop cuts org; %g; %g \n",feopMim,feopMax);
	163	//printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
	164
e156c3bb	165	Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
c2690925	166	AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
8c1c76e9	167	AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax));
c2690925	168	Int_t pdg = 0;
8c1c76e9	169	//if(eop>feopMim && eop<feopMax){
	170	//if(eop>feopMim && eop<feopMax){
	171	if(eop>feopMimCut && eop<feopMaxCut){
c2690925	172	pdg = 11;
8c1c76e9	173	//printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop);
e156c3bb	174	//if(pidqa)
e156c3bb	175	pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
c2690925	176	}
	177	else
	178	{
	179	pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
	180	}
	181
e156c3bb	182	AliDebug(1, Form("eID %g ; %d \n",eop,pdg));
8c1c76e9	183	//printf("eID %g ; %d \n",eop,pdg);
c2690925	184
	185	return pdg;
	186	}
e156c3bb	187
e156c3bb	188
8c1c76e9	189	//__________________________________________________________________________
	190	Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const
	191	{
	192	double maxCut = 0.0;
	193	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	194	if(esdtrack==NULL)return maxCut;
	195	double Pt = esdtrack->Pt();
	196
	197	if(flageop<0.5)
	198	{
	199	double meanP[3] = {0.991,1.0819,0.235};
	200	double sigP[3] = {6.52e-02,2.04e-01,3.34e-01};
	201	double mean = meanP[0]tanh(meanP[1]+meanP[2]Pt);
	202	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
	203	maxCut = mean+3.0*sig;
	204	}
	205	else
	206	{
	207	double meanP[3] = {0.99,1.299,0.35};
	208	double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
	209	double mean = meanP[0]tanh(meanP[1]+meanP[2]Pt);
	210	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
	211	maxCut = mean+3.0*sig;
	212	}
	213	return maxCut;
	214	}
	215
	216	Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
	217	{
	218	double mimCut = 0.0;
	219	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	220	if(esdtrack==NULL)return mimCut;
	221	double Pt = esdtrack->Pt();
	222
	223	if(flageop<0.5) // real
	224	{
	225	//printf("real");
	226	double meanP[3] = {0.991,1.0819,0.235};
	227	double sigP[3] = {6.52e-02,2.04e-01,3.34e-01};
	228	double mean = meanP[0]tanh(meanP[1]+meanP[2]Pt);
	229	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
	230	mimCut = mean-2.0*sig;
	231	}
	232	else // MC
	233	{
	234	//printf("MC");
	235	double meanP[3] = {0.99,1.299,0.35};
	236	double sigP[3] = {4.11e-02,1.588e-01,2.664e-01};
	237	double mean = meanP[0]tanh(meanP[1]+meanP[2]Pt);
	238	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*Pt);
	239	mimCut = mean-2.0*sig;
	240	}
	241	return mimCut;
	242	}
	243
e156c3bb	244	//___________________________________________________________________________
	245	Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
	246	{ // Returns e/p & Rmatch
	247
	248	Double_t matchclsE = 9999.9;
	249	double feop = -9999.9;
	250
	251	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
63bffcf1	252	if(esdtrack==NULL)return feop;
0ecbfc1b	253	const AliESDEvent *evt = esdtrack->GetESDEvent();
e156c3bb	254
8c1c76e9	255	Int_t icl = esdtrack->GetEMCALcluster();
e156c3bb	256
	257	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
	258	if(!cluster->IsEMCAL()) {return feop;}
	259
	260	matchclsE = cluster->E();
	261
	262
	263	if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
	264
	265	return feop;
	266
	267	}
	268
	269
c2690925	270	/*
e156c3bb	271	//____________________________________________________________________________________
e156c3bb	272	Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
c2690925	273	{ // Returns e/p if an electron is matched
	274
	275	Float_t clsPos[3];
	276	Double_t trkPos[3];
	277	Double_t matchclsE = 9999.9;
	278
	279	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	280	AliESDEvent *evt = esdtrack->GetESDEvent();
	281	Double_t magF = evt->GetMagneticField();
	282	Double_t magSign = 1.0;
	283	if(magF<0)magSign = -1.0;
	284	//printf("magF ; %g ; %g \n", magF,magSign);
	285
	286	if (!TGeoGlobalMagField::Instance()->GetField()) {
	287	printf("Loading field map...\n");
	288	//AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
	289	AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
	290	TGeoGlobalMagField::Instance()->SetField(field);
	291	}
	292
	293	AliEMCALTrack emctrack = new AliEMCALTrack(esdtrack);
	294	Double_t fieldB[3];
	295	emctrack->GetBxByBz(fieldB);
	296	//printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
	297
	298	for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
	299	{
	300
	301	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
	302	if(!cluster->IsEMCAL()) continue;
	303
	304	cluster->GetPosition(clsPos);
	305	if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
	306	emctrack->GetXYZ(trkPos);
	307
	308	TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
	309	TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
	310
	311	Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
	312	Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
	313
	314	double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
	315
	316	if(rmatch<0.02)
	317	{
	318	matchclsE = cluster->E();
	319	}
	320	}
	321	delete emctrack;
	322
	323	double feop = -9999.9;
	324	if(matchclsE<9999) feop = matchclsE/esdtrack->P();
	325
	326	// if(feop!=-9999.9)printf("%f\n",feop) ;
	327	return feop;
	328
	329	}
	330	*/
	331