[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEemcalPIDqa.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 AliHFEemcalPIDqa
//
// Monitoring EMCAL PID in the HFE PID montioring framework. The following
// quantities are monitored:
//   TPC Signal distribution 
// (Always as function of momentum, particle species and centrality 
// before and after cut)
// More information about the PID monitoring framework can be found in
// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
//
// Author:
//
//   Shingo Sakai <ssakai@lbl.gov>
#include <TClass.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>

#include <TMath.h>
#include "AliESDInputHandler.h"
//#include "AliVCluster.h"
//#include "AliVCaloCells.h"
//#include "AliVEvent.h"
#include "AliMagF.h"

#include "AliLog.h"
#include "AliPID.h"
#include "AliVParticle.h"
//#include "AliVTrack.h"
//#include "AliESDtrack.h"
#include "AliHFEcollection.h"
#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidEMCAL.h"
#include "AliHFEtools.h"
#include "AliHFEemcalPIDqa.h"

ClassImp(AliHFEemcalPIDqa)

//_________________________________________________________
AliHFEemcalPIDqa::AliHFEemcalPIDqa():
    AliHFEdetPIDqa()
  , fHistos(NULL)
{
  //
  // Dummy constructor
  //
}

//_________________________________________________________
AliHFEemcalPIDqa::AliHFEemcalPIDqa(const char* name):
    AliHFEdetPIDqa(name, "QA for EMCAL")
  , fHistos(NULL)
{
  //
  // Default constructor
  //
}

//_________________________________________________________
AliHFEemcalPIDqa::AliHFEemcalPIDqa(const AliHFEemcalPIDqa &o):
    AliHFEdetPIDqa(o)
  , fHistos()
{
  //
  // Copy constructor
  //
  o.Copy(*this);
}

//_________________________________________________________
AliHFEemcalPIDqa &AliHFEemcalPIDqa::operator=(const AliHFEemcalPIDqa &o){
  //
  // Do assignment
  //
  AliHFEdetPIDqa::operator=(o);
  if(&o != this) o.Copy(*this);
  return *this;
}

//_________________________________________________________
AliHFEemcalPIDqa::~AliHFEemcalPIDqa(){
  //
  // Destructor
  //
  if(fHistos) delete fHistos;
}

//_________________________________________________________
void AliHFEemcalPIDqa::Copy(TObject &o) const {
  //
  // Make copy
  //
  AliHFEemcalPIDqa &target = dynamic_cast<AliHFEemcalPIDqa &>(o);
  if(target.fHistos){
    delete target.fHistos;
    target.fHistos = NULL;
  }
  if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
}

//_________________________________________________________
Long64_t AliHFEemcalPIDqa::Merge(TCollection *coll){
  //
  // Merge with other objects
  //
  if(!coll) return 0;
  if(coll->IsEmpty()) return 1;

  TIter it(coll);
  AliHFEemcalPIDqa *refQA = NULL;
  TObject *o = NULL;
  Long64_t count = 0;
  TList listHistos;
  while((o = it())){
    refQA = dynamic_cast<AliHFEemcalPIDqa *>(o);
    if(!refQA) continue;

    listHistos.Add(refQA->fHistos);
    count++; 
  }
  fHistos->Merge(&listHistos);
  return count + 1;
}

//_________________________________________________________
void AliHFEemcalPIDqa::Initialize(){
  //
  // Define Histograms
  //

  fHistos = new AliHFEcollection("emcalqahistos", "Collection of EMCAL QA histograms");

  // Make common binning
  const Int_t kCentralityBins = 11;
  const Double_t kMinP = 0.;
  const Double_t kMaxP = 20.;
  const Double_t kTPCSigMim = 40.;
  const Double_t kTPCSigMax = 140.;

  // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, Centrality)
  Int_t nBins[6] = {AliPID::kSPECIES + 1, 20, 20, 400, kCentralityBins, 2};
  Double_t min[6] = {-1, kMinP, kMinP, kTPCSigMim,  0, 0.};
  Double_t max[6] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax,  11., 2.};
  fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; Centrality; PID Step", 6, nBins, min, max);

  //2nd histogram: EMCAL signal - E/p and Rmatch
  Int_t nBins2[6] = {AliPID::kSPECIES + 1, 40, 40, 1000, 250, 2};
  Double_t min2[6] = {-1, kMinP, kMinP, 0,  0, 0.};
  Double_t max2[6] = {AliPID::kSPECIES, kMaxP, kMaxP, 10,  1., 2.};
  fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; PID Step", 6, nBins2, min2, max2);
    
}


//_________________________________________________________
void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
  //
  // Fill TPC histograms
  //
  //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
  Float_t centrality = track->GetCentrality();

  //const AliVTrack *vtrack = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
  //const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(vtrack);
  const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track->GetRecTrack());

  Int_t species = track->GetAbInitioPID();
  if(species >= AliPID::kSPECIES) species = -1;

  Double_t contentSignal[6];
  contentSignal[0] = species;
  contentSignal[1] = track->GetRecTrack()->P();
  contentSignal[2] = track->GetRecTrack()->Pt();
  contentSignal[3] = esdtrack->GetTPCsignal(); //?
  contentSignal[4] = centrality;
  contentSignal[5] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;

  TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack);


  Double_t contentSignal2[6];
  contentSignal2[0] = species;
  contentSignal2[1] = track->GetRecTrack()->P();
  contentSignal2[2] = track->GetRecTrack()->Pt();
  contentSignal2[3] = emcsignal(0);
  contentSignal2[4] = emcsignal(1);
  contentSignal2[5] = contentSignal[5];

  //printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]); 
  fHistos->Fill("EMCAL_TPCdedx", contentSignal);
  fHistos->Fill("EMCAL_Signal", contentSignal2);
}

//_________________________________________________________
TH1 *AliHFEemcalPIDqa::GetHistogram(const char *name){
  // 
  // Get the histogram
  //
  if(!fHistos) return NULL;
  TObject *histo = fHistos->Get(name);
  if(!histo->InheritsFrom("TH1")) return NULL;
  return dynamic_cast<TH1 *>(histo);
}


//___________________________________________________________________________
TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const
{ // Returns e/p & Rmatch

  Float_t  clsPos[3];
  Double_t trkPos[3];
  Double_t matchclsE = 9999.9;
  TVector3 refVec(-9999,-9999,-9999);

  AliESDEvent *evt = esdtrack->GetESDEvent();

   //Int_t icl = esdtrack->GetEMCALcluster();
   Int_t icl = (const_cast<AliESDtrack *>(esdtrack))->GetEMCALcluster();

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

   cluster->GetPosition(clsPos);
   esdtrack->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));

   matchclsE = cluster->E();

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

   //   if(feop!=-9999.9)printf("%f\n",feop) ; 
   TVector3 emcsignal(feop,rmatch,0);

   return emcsignal;

}


/*
//___________________________________________________________________________
TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const
{ // Returns e/p & Rmatch

  Float_t  clsPos[3];
  Double_t trkPos[3];
  Double_t Rmatch;
  Double_t matchclsE = 9999.9;
  TVector3 refVec(-9999,-9999,-9999);

  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()) return refVec;

   cluster->GetPosition(clsPos);
   if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  return refVec;
   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();
       Rmatch = rmatch;
      }
  }
   delete emctrack;

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

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

}
*/
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 AliHFEemcalPIDqa
	17	//
	18	// Monitoring EMCAL PID in the HFE PID montioring framework. The following
	19	// quantities are monitored:
	20	// TPC Signal distribution
	21	// (Always as function of momentum, particle species and centrality
	22	// before and after cut)
	23	// More information about the PID monitoring framework can be found in
	24	// AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
	25	//
	26	// Author:
	27	//
	28	// Shingo Sakai <ssakai@lbl.gov>
	29	#include <TClass.h>
	30	#include <TH2.h>
	31	#include <THnSparse.h>
	32	#include <TString.h>
	33
e156c3bb	34	#include <TMath.h>
	35	#include "AliESDInputHandler.h"
	36	//#include "AliVCluster.h"
	37	//#include "AliVCaloCells.h"
	38	//#include "AliVEvent.h"
	39	#include "AliMagF.h"
	40
c2690925	41	#include "AliLog.h"
	42	#include "AliPID.h"
	43	#include "AliVParticle.h"
e156c3bb	44	//#include "AliVTrack.h"
e156c3bb	45	//#include "AliESDtrack.h"
c2690925	46	#include "AliHFEcollection.h"
	47	#include "AliHFEpid.h"
	48	#include "AliHFEpidBase.h"
	49	#include "AliHFEpidQAmanager.h"
	50	#include "AliHFEpidTPC.h"
	51	#include "AliHFEpidEMCAL.h"
	52	#include "AliHFEtools.h"
	53	#include "AliHFEemcalPIDqa.h"
	54
e156c3bb	55	ClassImp(AliHFEemcalPIDqa)
c2690925	56
	57	//_________________________________________________________
	58	AliHFEemcalPIDqa::AliHFEemcalPIDqa():
	59	AliHFEdetPIDqa()
	60	, fHistos(NULL)
	61	{
	62	//
	63	// Dummy constructor
	64	//
	65	}
	66
	67	//_________________________________________________________
	68	AliHFEemcalPIDqa::AliHFEemcalPIDqa(const char* name):
	69	AliHFEdetPIDqa(name, "QA for EMCAL")
	70	, fHistos(NULL)
	71	{
	72	//
	73	// Default constructor
	74	//
	75	}
	76
	77	//_________________________________________________________
	78	AliHFEemcalPIDqa::AliHFEemcalPIDqa(const AliHFEemcalPIDqa &o):
	79	AliHFEdetPIDqa(o)
	80	, fHistos()
	81	{
	82	//
	83	// Copy constructor
	84	//
	85	o.Copy(*this);
	86	}
	87
	88	//_________________________________________________________
	89	AliHFEemcalPIDqa &AliHFEemcalPIDqa::operator=(const AliHFEemcalPIDqa &o){
	90	//
	91	// Do assignment
	92	//
	93	AliHFEdetPIDqa::operator=(o);
	94	if(&o != this) o.Copy(*this);
	95	return *this;
	96	}
	97
	98	//_________________________________________________________
	99	AliHFEemcalPIDqa::~AliHFEemcalPIDqa(){
	100	//
	101	// Destructor
	102	//
	103	if(fHistos) delete fHistos;
	104	}
	105
	106	//_________________________________________________________
	107	void AliHFEemcalPIDqa::Copy(TObject &o) const {
	108	//
	109	// Make copy
	110	//
	111	AliHFEemcalPIDqa &target = dynamic_cast<AliHFEemcalPIDqa &>(o);
	112	if(target.fHistos){
	113	delete target.fHistos;
	114	target.fHistos = NULL;
	115	}
	116	if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
	117	}
	118
	119	//_________________________________________________________
120	Long64_t AliHFEemcalPIDqa::Merge(TCollection *coll){
121	//
122	// Merge with other objects
123	//
124	if(!coll) return 0;
125	if(coll->IsEmpty()) return 1;
126
127	TIter it(coll);
128	AliHFEemcalPIDqa *refQA = NULL;
129	TObject *o = NULL;
130	Long64_t count = 0;
131	TList listHistos;
132	while((o = it())){
133	refQA = dynamic_cast<AliHFEemcalPIDqa *>(o);
134	if(!refQA) continue;
135
136	listHistos.Add(refQA->fHistos);
137	count++;
138	}
139	fHistos->Merge(&listHistos);
140	return count + 1;
141	}
142
143	//_________________________________________________________
144	void AliHFEemcalPIDqa::Initialize(){
145	//
146	// Define Histograms
147	//
148
149	fHistos = new AliHFEcollection("emcalqahistos", "Collection of EMCAL QA histograms");
150
151	// Make common binning
152	const Int_t kCentralityBins = 11;
153	const Double_t kMinP = 0.;
154	const Double_t kMaxP = 20.;
155	const Double_t kTPCSigMim = 40.;
156	const Double_t kTPCSigMax = 140.;
157
158	// 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, Centrality)
e156c3bb	159	Int_t nBins[6] = {AliPID::kSPECIES + 1, 20, 20, 400, kCentralityBins, 2};
	160	Double_t min[6] = {-1, kMinP, kMinP, kTPCSigMim, 0, 0.};
	161	Double_t max[6] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 11., 2.};
	162	fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; Centrality; PID Step", 6, nBins, min, max);
	163
	164	//2nd histogram: EMCAL signal - E/p and Rmatch
	165	Int_t nBins2[6] = {AliPID::kSPECIES + 1, 40, 40, 1000, 250, 2};
	166	Double_t min2[6] = {-1, kMinP, kMinP, 0, 0, 0.};
	167	Double_t max2[6] = {AliPID::kSPECIES, kMaxP, kMaxP, 10, 1., 2.};
	168	fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; PID Step", 6, nBins2, min2, max2);
	169
c2690925	170	}
c2690925	171
e156c3bb	172
	173
	174
c2690925	175	//_________________________________________________________
e156c3bb	176	void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
c2690925	177	//
	178	// Fill TPC histograms
	179	//
	180	//AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
	181	Float_t centrality = track->GetCentrality();
	182
e156c3bb	183	//const AliVTrack vtrack = dynamic_cast<const AliVTrack >(track->GetRecTrack());
e156c3bb	184	//const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(vtrack);
6c70d827	185	const AliESDtrack esdtrack = static_cast<const AliESDtrack >(track->GetRecTrack());
c2690925	186
e156c3bb	187	Int_t species = track->GetAbInitioPID();
	188	if(species >= AliPID::kSPECIES) species = -1;
	189
	190	Double_t contentSignal[6];
	191	contentSignal[0] = species;
	192	contentSignal[1] = track->GetRecTrack()->P();
	193	contentSignal[2] = track->GetRecTrack()->Pt();
	194	contentSignal[3] = esdtrack->GetTPCsignal(); //?
	195	contentSignal[4] = centrality;
	196	contentSignal[5] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
	197
	198	TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack);
	199
	200
	201	Double_t contentSignal2[6];
	202	contentSignal2[0] = species;
	203	contentSignal2[1] = track->GetRecTrack()->P();
	204	contentSignal2[2] = track->GetRecTrack()->Pt();
	205	contentSignal2[3] = emcsignal(0);
	206	contentSignal2[4] = emcsignal(1);
	207	contentSignal2[5] = contentSignal[5];
	208
c2690925	209	//printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]);
c2690925	210	fHistos->Fill("EMCAL_TPCdedx", contentSignal);
e156c3bb	211	fHistos->Fill("EMCAL_Signal", contentSignal2);
c2690925	212	}
	213
	214	//_________________________________________________________
	215	TH1 AliHFEemcalPIDqa::GetHistogram(const char name){
	216	//
	217	// Get the histogram
	218	//
	219	if(!fHistos) return NULL;
	220	TObject *histo = fHistos->Get(name);
	221	if(!histo->InheritsFrom("TH1")) return NULL;
	222	return dynamic_cast<TH1 *>(histo);
	223	}
	224
e156c3bb	225
	226	//___________________________________________________________________________
	227	TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const
	228	{ // Returns e/p & Rmatch
	229
	230	Float_t clsPos[3];
	231	Double_t trkPos[3];
	232	Double_t matchclsE = 9999.9;
	233	TVector3 refVec(-9999,-9999,-9999);
	234
	235	AliESDEvent *evt = esdtrack->GetESDEvent();
	236
	237	//Int_t icl = esdtrack->GetEMCALcluster();
	238	Int_t icl = (const_cast<AliESDtrack *>(esdtrack))->GetEMCALcluster();
	239
	240	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
	241	if(!cluster->IsEMCAL()) return refVec;
	242
	243	cluster->GetPosition(clsPos);
	244	esdtrack->GetXYZ(trkPos);
	245
	246	TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
	247	TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
	248
	249
	250	Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
	251	Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
	252
	253	double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
	254
	255	matchclsE = cluster->E();
	256
	257	//double feop = -9999.9;
	258	//if(matchclsE<9999)
	259	double feop = matchclsE/esdtrack->P();
	260
	261	// if(feop!=-9999.9)printf("%f\n",feop) ;
	262	TVector3 emcsignal(feop,rmatch,0);
	263
	264	return emcsignal;
	265
	266	}
	267
	268
	269
	270
	271	/*
	272	//___________________________________________________________________________
	273	TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const
	274	{ // Returns e/p & Rmatch
	275
	276	Float_t clsPos[3];
	277	Double_t trkPos[3];
	278	Double_t Rmatch;
	279	Double_t matchclsE = 9999.9;
	280	TVector3 refVec(-9999,-9999,-9999);
	281
	282	AliESDEvent *evt = esdtrack->GetESDEvent();
	283	Double_t magF = evt->GetMagneticField();
	284	Double_t magSign = 1.0;
	285	if(magF<0)magSign = -1.0;
	286	//printf("magF ; %g ; %g \n", magF,magSign);
	287
	288	if (!TGeoGlobalMagField::Instance()->GetField()) {
289	printf("Loading field map...\n");
290	//AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
291	AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
292	TGeoGlobalMagField::Instance()->SetField(field);
293	}
294
295	AliEMCALTrack emctrack = new AliEMCALTrack(esdtrack);
296	Double_t fieldB[3];
297	emctrack->GetBxByBz(fieldB);
298	//printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
299
300	for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
301	{
302
303	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
304	if(!cluster->IsEMCAL()) return refVec;
305
306	cluster->GetPosition(clsPos);
307	if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) return refVec;
308	emctrack->GetXYZ(trkPos);
309
310	TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
311	TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
312
313	Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
314	Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
315
316	double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
317
318	if(rmatch<0.02)
319	{
320	matchclsE = cluster->E();
321	Rmatch = rmatch;
322	}
323	}
324	delete emctrack;
325
326	//double feop = -9999.9;
327	//if(matchclsE<9999)
328	double feop = matchclsE/esdtrack->P();
329
330	// if(feop!=-9999.9)printf("%f\n",feop) ;
331	TVector3 emcsignal(feop, Rmatch, 0);
332
333	return emcsignal;
334
335	}
336	*/
337