[u/mrichter/AliRoot.git] / PWGHF / 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.9)
  , feopMax(1.3)
{
  //
  // Constructor
  //
} 

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

//___________________________________________________________________
AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
  AliHFEpidBase("")
  , fPID(NULL)
  , feopMim(0.9)
  , feopMax(1.3)
{  
  // 
  // 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
{
 // max nSig cuts for 10d

  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)
    { //<--- new para for updated non-liniarity
     double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02}; 
     double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; 
     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] = {9.94183e-01,1.23952e+00,3.55571e-01}; 
     double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-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; 
    }

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

     return maxCut;
}

Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
{
 // mim nsig cuts for 10d
  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] = {1.08833e+00,1.17837e+00,7.47923e-02}; 
     double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; 
     double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
     double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
     mimCut = mean-3.0*sig;
    }
   else // MC
    { 
      //printf("MC"); 
     double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01}; 
     double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01}; 
     double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
     double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
     mimCut = mean-3.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==NULL) {return feop;}
   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)
fb87d707	50	, feopMim(0.9)
fb87d707	51	, feopMax(1.3)
c2690925	52	{
	53	//
	54	// Constructor
	55	//
	56	}
	57
	58	//___________________________________________________________________
	59	AliHFEpidEMCAL::AliHFEpidEMCAL(const Char_t *name):
	60	AliHFEpidBase(name)
	61	, fPID(NULL)
fb87d707	62	, feopMim(0.9)
fb87d707	63	, feopMax(1.3)
c2690925	64	{
	65	//
	66	// Constructor
	67	//
	68	}
	69
	70	//___________________________________________________________________
	71	AliHFEpidEMCAL::AliHFEpidEMCAL(const AliHFEpidEMCAL &c):
	72	AliHFEpidBase("")
	73	, fPID(NULL)
fb87d707	74	, feopMim(0.9)
fb87d707	75	, feopMax(1.3)
c2690925	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	//
fb87d707	128
63bffcf1	129	if(track==NULL)return 0;
63bffcf1	130
8c1c76e9	131	if(!fkPIDResponse) return 0;
c2690925	132	AliDebug(2, "PID object available");
c2690925	133	// EMCal not fESDpid (s.s Feb. 11)
63bffcf1	134
	135	const AliVTrack trk = dynamic_cast<const AliVTrack >(track->GetRecTrack());
	136	if (trk == NULL) return 0;
e156c3bb	137	//AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
63bffcf1	138	if(!(trk->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
c2690925	139	AliDebug(2, "Track Has EMCAL PID");
e156c3bb	140
c9f1d8ef	141	if(pidqa)
e156c3bb	142	pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
c2690925	143	// not QA for EMCal, will be added (s.s Feb. 11)
c2690925	144
8c1c76e9	145
	146	double feopMimCut = 0.0;
	147	double feopMaxCut = 0.0;
	148
	149	feopMimCut = feopMim;
	150	feopMaxCut = feopMax;
	151
	152	if(feopMax >900) // not use fix cut
	153	{
	154	feopMimCut = CalEopCutMim(track->GetRecTrack(),0);
	155	feopMaxCut = CalEopCutMax(track->GetRecTrack(),0);
	156	}
	157	if(feopMax < -900) // not use fix cut for MC
	158	{
	159	feopMimCut = CalEopCutMim(track->GetRecTrack(),1);
	160	feopMaxCut = CalEopCutMax(track->GetRecTrack(),1);
	161	}
	162
a8ef1999	163	//printf("eop cuts org; %g; %g \n",feopMim,feopMax);
a8ef1999	164	//printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
8c1c76e9	165
e156c3bb	166	Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
c2690925	167	AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
8c1c76e9	168	AliDebug(2, Form("E/p cut ; %g ; %g \n", feopMim,feopMax));
c2690925	169	Int_t pdg = 0;
8c1c76e9	170	//if(eop>feopMim && eop<feopMax){
	171	//if(eop>feopMim && eop<feopMax){
	172	if(eop>feopMimCut && eop<feopMaxCut){
c2690925	173	pdg = 11;
8c1c76e9	174	//printf("eop cuts ; %g; %g ; %g\n",feopMimCut,feopMaxCut,eop);
c9f1d8ef	175	if(pidqa)
e156c3bb	176	pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
c2690925	177	}
	178	else
	179	{
	180	pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
	181	}
	182
e156c3bb	183	AliDebug(1, Form("eID %g ; %d \n",eop,pdg));
8c1c76e9	184	//printf("eID %g ; %d \n",eop,pdg);
c2690925	185
	186	return pdg;
	187	}
e156c3bb	188
e156c3bb	189
8c1c76e9	190	//__________________________________________________________________________
	191	Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t flageop) const
	192	{
57189f04	193	// max nSig cuts for 10d
57189f04	194
8c1c76e9	195	double maxCut = 0.0;
	196	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	197	if(esdtrack==NULL)return maxCut;
57189f04	198	double pt = esdtrack->Pt();
8c1c76e9	199
8c1c76e9	200	if(flageop<0.5)
c9f1d8ef	201	{ //<--- new para for updated non-liniarity
a8ef1999	202	double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02};
a8ef1999	203	double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04};
57189f04	204	double mean = meanP[0]tanh(meanP[1]+meanP[2]pt);
57189f04	205	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
8c1c76e9	206	maxCut = mean+3.0*sig;
	207	}
	208	else
	209	{
a8ef1999	210	double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01};
a8ef1999	211	double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01};
57189f04	212	double mean = meanP[0]tanh(meanP[1]+meanP[2]pt);
57189f04	213	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
8c1c76e9	214	maxCut = mean+3.0*sig;
8c1c76e9	215	}
a8ef1999	216
	217	//printf("eop cuts org; %g; %g ; %g\n",feopMim,feopMax,pt);
	218
8c1c76e9	219	return maxCut;
	220	}
	221
	222	Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t flageop) const
	223	{
57189f04	224	// mim nsig cuts for 10d
8c1c76e9	225	double mimCut = 0.0;
	226	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	227	if(esdtrack==NULL)return mimCut;
57189f04	228	double pt = esdtrack->Pt();
8c1c76e9	229
	230	if(flageop<0.5) // real
	231	{
	232	//printf("real");
a8ef1999	233	double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02};
a8ef1999	234	double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04};
57189f04	235	double mean = meanP[0]tanh(meanP[1]+meanP[2]pt);
57189f04	236	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
a8ef1999	237	mimCut = mean-3.0*sig;
8c1c76e9	238	}
	239	else // MC
	240	{
	241	//printf("MC");
a8ef1999	242	double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01};
a8ef1999	243	double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01};
57189f04	244	double mean = meanP[0]tanh(meanP[1]+meanP[2]pt);
57189f04	245	double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt);
a8ef1999	246	mimCut = mean-3.0*sig;
8c1c76e9	247	}
	248	return mimCut;
	249	}
	250
e156c3bb	251	//___________________________________________________________________________
	252	Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
	253	{ // Returns e/p & Rmatch
	254
	255	Double_t matchclsE = 9999.9;
	256	double feop = -9999.9;
	257
	258	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
63bffcf1	259	if(esdtrack==NULL)return feop;
0ecbfc1b	260	const AliESDEvent *evt = esdtrack->GetESDEvent();
e156c3bb	261
8c1c76e9	262	Int_t icl = esdtrack->GetEMCALcluster();
e156c3bb	263
e156c3bb	264	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
3a7690e5	265	if(cluster==NULL) {return feop;}
e156c3bb	266	if(!cluster->IsEMCAL()) {return feop;}
	267
	268	matchclsE = cluster->E();
e156c3bb	269
	270	if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
	271
	272	return feop;
	273
	274	}
	275
	276
c2690925	277	/*
e156c3bb	278	//____________________________________________________________________________________
e156c3bb	279	Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
c2690925	280	{ // Returns e/p if an electron is matched
	281
	282	Float_t clsPos[3];
	283	Double_t trkPos[3];
	284	Double_t matchclsE = 9999.9;
	285
	286	const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(track);
	287	AliESDEvent *evt = esdtrack->GetESDEvent();
	288	Double_t magF = evt->GetMagneticField();
	289	Double_t magSign = 1.0;
	290	if(magF<0)magSign = -1.0;
	291	//printf("magF ; %g ; %g \n", magF,magSign);
	292
	293	if (!TGeoGlobalMagField::Instance()->GetField()) {
	294	printf("Loading field map...\n");
	295	//AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
	296	AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
	297	TGeoGlobalMagField::Instance()->SetField(field);
	298	}
	299
	300	AliEMCALTrack emctrack = new AliEMCALTrack(esdtrack);
	301	Double_t fieldB[3];
	302	emctrack->GetBxByBz(fieldB);
	303	//printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
	304
	305	for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
	306	{
	307
	308	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
	309	if(!cluster->IsEMCAL()) continue;
	310
	311	cluster->GetPosition(clsPos);
	312	if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
	313	emctrack->GetXYZ(trkPos);
	314
	315	TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
	316	TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
	317
	318	Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
	319	Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
	320
	321	double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
	322
	323	if(rmatch<0.02)
	324	{
	325	matchclsE = cluster->E();
	326	}
	327	}
	328	delete emctrack;
	329
	330	double feop = -9999.9;
	331	if(matchclsE<9999) feop = matchclsE/esdtrack->P();
	332
	333	// if(feop!=-9999.9)printf("%f\n",feop) ;
	334	return feop;
	335
	336	}
	337	*/
	338