[u/mrichter/AliRoot.git] / PWGHF / 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 "AliPIDResponse.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"
#include "AliTrackerBase.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 = 50.;
  const Double_t kTPCSigMim = 40.;
  const Double_t kTPCSigMax = 140.;

  // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, Centrality, select)
  Int_t nBins[12] = {AliPID::kSPECIES + 1, 500, 500,          400, 630,   200,   600,  300, 100, kCentralityBins, 6, 2};
  Double_t min[12] = {-1,               kMinP, kMinP,  kTPCSigMim,  0.,  -1.0,  -8.0,    0,   0,            0, -3.0, 0.};
  Double_t max[12] = {AliPID::kSPECIES, kMaxP, kMaxP,  kTPCSigMax, 6.3,   1.0,   4.0,  3.0, 0.1,           11., 3.0, 2.};
  fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ;Centrality; PID Step; ", 12, nBins, min, max);

  //2nd histogram: EMCAL signal - E/p 
  Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 100, 600, 2};
  Double_t min2[7] = {-1, kMinP, kMinP, 0,  0, -8., 0.};
  Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5,  0.1, 4., 2.};
  fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; TPCnsigma; PID Step", 7, 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;

  // Get TPC nSigma (based on AliHFEtpcPIDqa)
  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid)); 
  const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
  Float_t nSigmatpc = pidResponse ? pidResponse->NumberOfSigmasTPC(static_cast<const AliVTrack *>(track->GetRecTrack()), AliPID::kElectron) : 0.; 
  //printf("nSigmatpc = %f\n",nSigmatpc);
  AliDebug(2, Form("nSigmatpc = %f\n",nSigmatpc));
  //

  double charge = esdtrack->Charge();
  //printf("charge %f\n",charge); 


  //printf("phi %f, eta %f\n;",phi,eta);

  // Get E/p
  TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack);

  Double_t contentSignal2[7];
  contentSignal2[0] = species;
  contentSignal2[1] = track->GetRecTrack()->P();
  contentSignal2[2] = track->GetRecTrack()->Pt();
  contentSignal2[3] = emcsignal(0);//e over p
  contentSignal2[4] = emcsignal(1);//residual
  contentSignal2[5] = nSigmatpc;
  contentSignal2[6] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;

  // QA array
  Double_t contentSignal[12];
  contentSignal[0] = species;
  contentSignal[1] = track->GetRecTrack()->P();
  contentSignal[2] = track->GetRecTrack()->Pt();
  contentSignal[3] = esdtrack->GetTPCsignal(); //?
  double phi  = track->GetRecTrack()->Phi();
  double eta  = track->GetRecTrack()->Eta();
  contentSignal[4] = phi;
  contentSignal[5] = eta;
  contentSignal[6] = nSigmatpc;
  contentSignal[7] = emcsignal(0);
  contentSignal[8] = emcsignal(2);
  contentSignal[9] = centrality;
  contentSignal[10] = charge;
  contentSignal[11] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;


  //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;
  Double_t fMass=0.139;
  Double_t fStep=10; //This is taken from EMCAL tender, hardcoded!
  TVector3 refVec(-9999,-9999,-9999);

  const AliESDEvent *evt = esdtrack->GetESDEvent();

  //Proper trackmatching, added by Tomas
  //Get matched cluster
  Int_t icl = esdtrack->GetEMCALcluster(); //From tender
  AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
  if(!cluster->IsEMCAL()) return refVec;

  // from ESDs
  double delphi = cluster->GetTrackDx(); 
  double deleta = cluster->GetTrackDz(); 
  double rmatch1 = sqrt(pow(delphi,2)+pow(deleta,2));

  //
  cluster->GetPosition(clsPos);
  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);//Vector of emcal cluster

  //Extrapolate track to emcal, properly! First, get external params!
  AliExternalTrackParam *trackParam = 0;
  const AliESDfriendTrack*  friendTrack = esdtrack->GetFriendTrack();
  if(friendTrack && friendTrack->GetTPCOut())
    {
      //Use TPC Out as starting point if it is available
      trackParam=  const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
    }
  else
    {
      //Otherwise use TPC inner
      trackParam =  const_cast<AliExternalTrackParam*>(esdtrack->GetInnerParam());
    }

  Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
  vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
  trackParam->Rotate(alpha); //Rotate the track to the same local extrapolation system 

  //Do extrapolation:
  if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return refVec; 
  trackParam->GetXYZ(trkPos); //Get the extrapolated global position, done!


  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,rmatch1);
  
  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"
8c1c76e9	40	#include "AliPIDResponse.h"
e156c3bb	41
c2690925	42	#include "AliLog.h"
	43	#include "AliPID.h"
	44	#include "AliVParticle.h"
e156c3bb	45	//#include "AliVTrack.h"
e156c3bb	46	//#include "AliESDtrack.h"
c2690925	47	#include "AliHFEcollection.h"
	48	#include "AliHFEpid.h"
	49	#include "AliHFEpidBase.h"
	50	#include "AliHFEpidQAmanager.h"
	51	#include "AliHFEpidTPC.h"
	52	#include "AliHFEpidEMCAL.h"
	53	#include "AliHFEtools.h"
	54	#include "AliHFEemcalPIDqa.h"
8c1c76e9	55	#include "AliTrackerBase.h"
c2690925	56
e156c3bb	57	ClassImp(AliHFEemcalPIDqa)
c2690925	58
	59	//_________________________________________________________
	60	AliHFEemcalPIDqa::AliHFEemcalPIDqa():
	61	AliHFEdetPIDqa()
	62	, fHistos(NULL)
	63	{
	64	//
	65	// Dummy constructor
	66	//
	67	}
	68
	69	//_________________________________________________________
	70	AliHFEemcalPIDqa::AliHFEemcalPIDqa(const char* name):
	71	AliHFEdetPIDqa(name, "QA for EMCAL")
	72	, fHistos(NULL)
	73	{
	74	//
	75	// Default constructor
	76	//
	77	}
	78
	79	//_________________________________________________________
	80	AliHFEemcalPIDqa::AliHFEemcalPIDqa(const AliHFEemcalPIDqa &o):
	81	AliHFEdetPIDqa(o)
	82	, fHistos()
	83	{
	84	//
	85	// Copy constructor
	86	//
	87	o.Copy(*this);
	88	}
	89
	90	//_________________________________________________________
	91	AliHFEemcalPIDqa &AliHFEemcalPIDqa::operator=(const AliHFEemcalPIDqa &o){
	92	//
	93	// Do assignment
	94	//
	95	AliHFEdetPIDqa::operator=(o);
	96	if(&o != this) o.Copy(*this);
	97	return *this;
	98	}
	99
	100	//_________________________________________________________
	101	AliHFEemcalPIDqa::~AliHFEemcalPIDqa(){
	102	//
	103	// Destructor
	104	//
	105	if(fHistos) delete fHistos;
	106	}
	107
	108	//_________________________________________________________
	109	void AliHFEemcalPIDqa::Copy(TObject &o) const {
	110	//
	111	// Make copy
	112	//
	113	AliHFEemcalPIDqa &target = dynamic_cast<AliHFEemcalPIDqa &>(o);
	114	if(target.fHistos){
	115	delete target.fHistos;
	116	target.fHistos = NULL;
	117	}
	118	if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
	119	}
	120
	121	//_________________________________________________________
122	Long64_t AliHFEemcalPIDqa::Merge(TCollection *coll){
123	//
124	// Merge with other objects
125	//
126	if(!coll) return 0;
127	if(coll->IsEmpty()) return 1;
128
129	TIter it(coll);
130	AliHFEemcalPIDqa *refQA = NULL;
131	TObject *o = NULL;
132	Long64_t count = 0;
133	TList listHistos;
134	while((o = it())){
135	refQA = dynamic_cast<AliHFEemcalPIDqa *>(o);
136	if(!refQA) continue;
137
138	listHistos.Add(refQA->fHistos);
139	count++;
140	}
141	fHistos->Merge(&listHistos);
142	return count + 1;
143	}
144
145	//_________________________________________________________
146	void AliHFEemcalPIDqa::Initialize(){
147	//
148	// Define Histograms
149	//
150
151	fHistos = new AliHFEcollection("emcalqahistos", "Collection of EMCAL QA histograms");
152
153	// Make common binning
154	const Int_t kCentralityBins = 11;
155	const Double_t kMinP = 0.;
8c1c76e9	156	const Double_t kMaxP = 50.;
c2690925	157	const Double_t kTPCSigMim = 40.;
	158	const Double_t kTPCSigMax = 140.;
	159
e17c1f86	160	// 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, Centrality, select)
a8ef1999	161	Int_t nBins[12] = {AliPID::kSPECIES + 1, 500, 500, 400, 630, 200, 600, 300, 100, kCentralityBins, 6, 2};
e17c1f86	162	Double_t min[12] = {-1, kMinP, kMinP, kTPCSigMim, 0., -1.0, -8.0, 0, 0, 0, -3.0, 0.};
a8ef1999	163	Double_t max[12] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 6.3, 1.0, 4.0, 3.0, 0.1, 11., 3.0, 2.};
e17c1f86	164	fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ;Centrality; PID Step; ", 12, nBins, min, max);
e156c3bb	165
8c1c76e9	166	//2nd histogram: EMCAL signal - E/p
a8ef1999	167	Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 100, 600, 2};
e17c1f86	168	Double_t min2[7] = {-1, kMinP, kMinP, 0, 0, -8., 0.};
a8ef1999	169	Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5, 0.1, 4., 2.};
8c1c76e9	170	fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; TPCnsigma; PID Step", 7, nBins2, min2, max2);
e156c3bb	171
c2690925	172	}
c2690925	173
e156c3bb	174
	175
	176
c2690925	177	//_________________________________________________________
e156c3bb	178	void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
c2690925	179	//
	180	// Fill TPC histograms
	181	//
	182	//AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
	183	Float_t centrality = track->GetCentrality();
	184
e156c3bb	185	//const AliVTrack vtrack = dynamic_cast<const AliVTrack >(track->GetRecTrack());
e156c3bb	186	//const AliESDtrack esdtrack = dynamic_cast<const AliESDtrack >(vtrack);
6c70d827	187	const AliESDtrack esdtrack = static_cast<const AliESDtrack >(track->GetRecTrack());
c2690925	188
e156c3bb	189	Int_t species = track->GetAbInitioPID();
	190	if(species >= AliPID::kSPECIES) species = -1;
	191
8c1c76e9	192	// Get TPC nSigma (based on AliHFEtpcPIDqa)
	193	AliHFEpidTPC tpcpid = dynamic_cast<AliHFEpidTPC >(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
	194	const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
	195	Float_t nSigmatpc = pidResponse ? pidResponse->NumberOfSigmasTPC(static_cast<const AliVTrack *>(track->GetRecTrack()), AliPID::kElectron) : 0.;
	196	//printf("nSigmatpc = %f\n",nSigmatpc);
	197	AliDebug(2, Form("nSigmatpc = %f\n",nSigmatpc));
	198	//
	199
e17c1f86	200	double charge = esdtrack->Charge();
e17c1f86	201	//printf("charge %f\n",charge);
c9f1d8ef	202
	203
	204	//printf("phi %f, eta %f\n;",phi,eta);
e156c3bb	205
8c1c76e9	206	// Get E/p
e156c3bb	207	TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack);
e156c3bb	208
8c1c76e9	209	Double_t contentSignal2[7];
e156c3bb	210	contentSignal2[0] = species;
	211	contentSignal2[1] = track->GetRecTrack()->P();
	212	contentSignal2[2] = track->GetRecTrack()->Pt();
8c1c76e9	213	contentSignal2[3] = emcsignal(0);//e over p
	214	contentSignal2[4] = emcsignal(1);//residual
	215	contentSignal2[5] = nSigmatpc;
c9f1d8ef	216	contentSignal2[6] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
	217
	218	// QA array
e17c1f86	219	Double_t contentSignal[12];
c9f1d8ef	220	contentSignal[0] = species;
	221	contentSignal[1] = track->GetRecTrack()->P();
	222	contentSignal[2] = track->GetRecTrack()->Pt();
	223	contentSignal[3] = esdtrack->GetTPCsignal(); //?
	224	double phi = track->GetRecTrack()->Phi();
	225	double eta = track->GetRecTrack()->Eta();
	226	contentSignal[4] = phi;
	227	contentSignal[5] = eta;
	228	contentSignal[6] = nSigmatpc;
	229	contentSignal[7] = emcsignal(0);
e17c1f86	230	contentSignal[8] = emcsignal(2);
	231	contentSignal[9] = centrality;
	232	contentSignal[10] = charge;
	233	contentSignal[11] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
c9f1d8ef	234
e156c3bb	235
c2690925	236	//printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]);
c2690925	237	fHistos->Fill("EMCAL_TPCdedx", contentSignal);
e156c3bb	238	fHistos->Fill("EMCAL_Signal", contentSignal2);
c2690925	239	}
	240
	241	//_________________________________________________________
	242	TH1 AliHFEemcalPIDqa::GetHistogram(const char name){
	243	//
	244	// Get the histogram
	245	//
	246	if(!fHistos) return NULL;
	247	TObject *histo = fHistos->Get(name);
	248	if(!histo->InheritsFrom("TH1")) return NULL;
	249	return dynamic_cast<TH1 *>(histo);
	250	}
	251
e156c3bb	252
	253	//___________________________________________________________________________
	254	TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const
	255	{ // Returns e/p & Rmatch
	256
	257	Float_t clsPos[3];
	258	Double_t trkPos[3];
	259	Double_t matchclsE = 9999.9;
8c1c76e9	260	Double_t fMass=0.139;
8c1c76e9	261	Double_t fStep=10; //This is taken from EMCAL tender, hardcoded!
e156c3bb	262	TVector3 refVec(-9999,-9999,-9999);
e156c3bb	263
0ecbfc1b	264	const AliESDEvent *evt = esdtrack->GetESDEvent();
e156c3bb	265
8c1c76e9	266	//Proper trackmatching, added by Tomas
	267	//Get matched cluster
	268	Int_t icl = esdtrack->GetEMCALcluster(); //From tender
	269	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
	270	if(!cluster->IsEMCAL()) return refVec;
e17c1f86	271
	272	// from ESDs
	273	double delphi = cluster->GetTrackDx();
	274	double deleta = cluster->GetTrackDz();
	275	double rmatch1 = sqrt(pow(delphi,2)+pow(deleta,2));
	276
	277	//
8c1c76e9	278	cluster->GetPosition(clsPos);
	279	TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);//Vector of emcal cluster
	280
	281	//Extrapolate track to emcal, properly! First, get external params!
	282	AliExternalTrackParam *trackParam = 0;
	283	const AliESDfriendTrack* friendTrack = esdtrack->GetFriendTrack();
	284	if(friendTrack && friendTrack->GetTPCOut())
	285	{
	286	//Use TPC Out as starting point if it is available
	287	trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
	288	}
	289	else
	290	{
	291	//Otherwise use TPC inner
	292	trackParam = const_cast<AliExternalTrackParam*>(esdtrack->GetInnerParam());
	293	}
	294
	295	Double_t alpha = ((int)(vec.Phi()TMath::RadToDeg()/20)+0.5)20*TMath::DegToRad();
	296	vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
	297	trackParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
	298
	299	//Do extrapolation:
	300	if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return refVec;
	301	trackParam->GetXYZ(trkPos); //Get the extrapolated global position, done!
	302
	303
	304
	305	TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
	306	TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
	307
	308
	309	Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
	310	Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
	311
	312	double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
	313
e17c1f86	314
8c1c76e9	315	matchclsE = cluster->E();
	316
	317	//double feop = -9999.9;
	318	//if(matchclsE<9999)
	319	double feop = matchclsE/esdtrack->P();
	320
	321	// if(feop!=-9999.9)printf("%f\n",feop) ;
e17c1f86	322	TVector3 emcsignal(feop,rmatch,rmatch1);
8c1c76e9	323
	324	return emcsignal;
	325
e156c3bb	326	}
	327
	328
	329
	330
	331	/*
	332	//___________________________________________________________________________
	333	TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const
	334	{ // Returns e/p & Rmatch
	335
	336	Float_t clsPos[3];
	337	Double_t trkPos[3];
	338	Double_t Rmatch;
	339	Double_t matchclsE = 9999.9;
	340	TVector3 refVec(-9999,-9999,-9999);
	341
	342	AliESDEvent *evt = esdtrack->GetESDEvent();
	343	Double_t magF = evt->GetMagneticField();
	344	Double_t magSign = 1.0;
	345	if(magF<0)magSign = -1.0;
	346	//printf("magF ; %g ; %g \n", magF,magSign);
	347
	348	if (!TGeoGlobalMagField::Instance()->GetField()) {
	349	printf("Loading field map...\n");
	350	//AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
	351	AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
	352	TGeoGlobalMagField::Instance()->SetField(field);
	353	}
	354
	355	AliEMCALTrack emctrack = new AliEMCALTrack(esdtrack);
	356	Double_t fieldB[3];
	357	emctrack->GetBxByBz(fieldB);
	358	//printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
	359
	360	for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
	361	{
	362
	363	AliVCluster cluster = (AliVCluster) evt->GetCaloCluster(icl);
	364	if(!cluster->IsEMCAL()) return refVec;
	365
	366	cluster->GetPosition(clsPos);
	367	if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) return refVec;
	368	emctrack->GetXYZ(trkPos);
	369
	370	TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
	371	TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
	372
	373	Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
	374	Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
	375
	376	double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
	377
	378	if(rmatch<0.02)
	379	{
	380	matchclsE = cluster->E();
	381	Rmatch = rmatch;
	382	}
	383	}
	384	delete emctrack;
	385
	386	//double feop = -9999.9;
	387	//if(matchclsE<9999)
	388	double feop = matchclsE/esdtrack->P();
	389
390	// if(feop!=-9999.9)printf("%f\n",feop) ;
391	TVector3 emcsignal(feop, Rmatch, 0);
392
393	return emcsignal;
394
395	}
396	*/
397