[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisEtReconstructed.cxx

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Base class for ESD analysis
//  - reconstruction output
//  implementation file
//
//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
//_________________________________________________________________________

#include "AliAnalysisEtReconstructed.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
#include "AliEMCALTrack.h"
#include "AliESDCaloCluster.h"
#include "TVector3.h"
#include "TGeoGlobalMagField.h"
#include "AliMagF.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliESDtrackCuts.h"
#include "AliVParticle.h"
#include "TDatabasePDG.h"
#include "TList.h"
#include "AliESDpid.h"
#include <iostream>
#include "TH2F.h"

using namespace std;

ClassImp(AliAnalysisEtReconstructed);


AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
        AliAnalysisEt()
        ,fTrackDistanceCut(0)
        ,fPidCut(0)
        ,fClusterType(0)
        ,fHistChargedPionEnergyDeposit(0)
        ,fHistProtonEnergyDeposit(0)
        ,fHistAntiProtonEnergyDeposit(0)
        ,fHistChargedKaonEnergyDeposit(0)
        ,fHistMuonEnergyDeposit(0)
{

}

AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed() 
{
}

Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
{ // analyse ESD event
    ResetEventValues();
    AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);

    Double_t protonMass = fPdgDB->GetParticle("proton")->Mass();

    //for PID
    AliESDpid *pID = new AliESDpid();
    pID->MakePID(event);
    TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event);
    Int_t nGoodTracks = list->GetEntries();
    // printf("nGoodTracks %d nCaloClusters %d\n", nGoodTracks, event->GetNumberOfCaloClusters());

    for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
      {
	AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
        if (!track)
        {
            Printf("ERROR: Could not get track %d", iTrack);
            continue;
        }

        fMultiplicity++;


	Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
	nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
	nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
	nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
	nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
	/*
	bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
	bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
	bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
	bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
	*/

        Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
        Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);

        Float_t massPart = 0;

        const Double_t *pidWeights = track->PID();
	Int_t maxpid = -1;
        Double_t maxpidweight = 0;
            
        if (pidWeights)
        {
            for (Int_t p =0; p < AliPID::kSPECIES; p++)
            {
                if (pidWeights[p] > maxpidweight)
                {
                    maxpidweight = pidWeights[p];
                    maxpid = p;
                }
            }
            if (maxpid == AliPID::kProton)
            {
	      //by definition of ET
		massPart = -protonMass*track->Charge();
            }

        }

        Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
	//printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout

        if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
        {
	    fTotChargedEt +=  et;
            fChargedMultiplicity++;
	    if (maxpid != -1)
            {
                if (maxpid == AliPID::kPion)
                {
                    fProtonEt += et;
                }
                if (maxpid == AliPID::kKaon)
                {
                    fChargedKaonEt += et;
                }
                if (maxpid == AliPID::kMuon)
                {
                    fMuonEt += et;
                }
                if (maxpid == AliPID::kElectron)
                {
                    fElectronEt += et;
                }
            }

            if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
            {
                fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
		if (maxpid != -1)
                {
                    if (maxpid == AliPID::kPion)
                    {
                        fProtonEtAcc += et;
                    }
                    if (maxpid == AliPID::kKaon)
                    {
                        fChargedKaonEtAcc += et;
                    }
                    if (maxpid == AliPID::kMuon)
                    {
                        fMuonEtAcc += et;
                    }
                    if (maxpid == AliPID::kElectron)
                    {
                        fElectronEtAcc += et;
                    }
                }
           
            }
        }

        if (TrackHitsCalorimeter(track, event->GetMagneticField()))
        {
	  Double_t phi = track->Phi();
	  Double_t pt = track->Pt();
	  // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
	  if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
	  else fHistPhivsPtNeg->Fill(phi, pt);
        }
      }
    delete pID;

    for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
    {
        AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
        if (!cluster)
        {
            Printf("ERROR: Could not get cluster %d", iCluster);
            continue;
        }

	if (cluster->GetType() != fClusterType) continue;
	//if(cluster->GetTracksMatched() > 0) 
	// printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout
	       
        
        if (cluster->E() < fClusterEnergyCut) continue;
        Float_t pos[3];
        TVector3 cp(pos);
        cluster->GetPosition(pos);

	Double_t distance = cluster->GetEmcCpvDistance();
	Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
	if ( cluster->IsEMCAL() ) {
	  distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
	}

	fHistTMDeltaR->Fill(distance);
        if (distance < fTrackDistanceCut)
        {
            if (cluster->GetNTracksMatched() == 1 && trackMatchedIndex>0)
            {
                AliVTrack *track = event->GetTrack(trackMatchedIndex);
                const Double_t *pidWeights = track->PID();
                
		Double_t maxpidweight = 0;
		Int_t maxpid = 0;
                
		if (pidWeights)
                {
                    for (Int_t p =0; p < AliPID::kSPECIES; p++)
                    {
                        if (pidWeights[p] > maxpidweight)
                        {
                            maxpidweight = pidWeights[p];
                            maxpid = p;
                        }
                    }
                    if(maxpidweight > fPidCut)
		    {
		       if(maxpid == AliPID::kProton)
		       {
			  if(track->Charge() == 1)
			  {
			     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
			  }
			  else if(track->Charge() == -1)
			  {
			     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
			  }
		       }
		       else if(maxpid == AliPID::kPion)
		       {
			  fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
		       }
		       else if(maxpid == AliPID::kKaon)
		       {
			  fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
		       }   
		       else if(maxpid == AliPID::kMuon)
		       {
			  fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
		       }
		    }
                }
            }

            continue;
        } // distance

        if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;

        cluster->GetPosition(pos);
      
	// TODO: replace with TVector3, too lazy now...

        float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);

        float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
        // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
        fTotNeutralEt += cluster->E() * TMath::Sin(theta);
	fNeutralMultiplicity++;

        fMultiplicity++;
    }

    fTotNeutralEtAcc = fTotNeutralEt;
    fTotEt = fTotChargedEt + fTotNeutralEt;
    fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;

    // Fill the histograms...
    FillHistograms();

    return 0;
}

bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
{ // check vertex

    Float_t bxy = 999.;
    Float_t bz = 999.;
    dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);

    bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) && 
      (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) && 
      (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) && 
      (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) && 
      (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut()); 

    return status;
}

void AliAnalysisEtReconstructed::Init()
{ // Init
    AliAnalysisEt::Init();
    fPidCut = fCuts->GetReconstructedPidCut();
    TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));

}

bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
{ // propagate track to detector radius

   AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
   // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());

    Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);

    // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
    return prop && 
		   TMath::Abs(esdTrack->Eta()) < fEtaCutAcc && 
		   esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. && 
		   esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
}

void AliAnalysisEtReconstructed::FillOutputList(TList* list)
{ // add some extra histograms to the ones from base class
    AliAnalysisEt::FillOutputList(list);

    list->Add(fHistChargedPionEnergyDeposit);
    list->Add(fHistProtonEnergyDeposit);
    list->Add(fHistAntiProtonEnergyDeposit);
    list->Add(fHistChargedKaonEnergyDeposit);
    list->Add(fHistMuonEnergyDeposit);
}

void AliAnalysisEtReconstructed::CreateHistograms()
{ // add some extra histograms to the ones from base class
    AliAnalysisEt::CreateHistograms();

    Int_t nbinsEt = 1000;
    Double_t minEt = 0;
    Double_t maxEt = 10;

    // possibly change histogram limits
    if (fCuts) {
      nbinsEt = fCuts->GetHistNbinsParticleEt();
      minEt = fCuts->GetHistMinParticleEt();
      maxEt = fCuts->GetHistMaxParticleEt();
    }

    TString histname;
    histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
    fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
    fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
    fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
    
    histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
    fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
    fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
    fHistProtonEnergyDeposit->SetYTitle("Energy of track");
    
    histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
    fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
    fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
    fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
    
    histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
    fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
    fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
    fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
    
    histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
    fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
    fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
    fHistMuonEnergyDeposit->SetYTitle("Energy of track");
    

}

Double_t 
AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
						     Int_t *trkMatchId,
						     const AliESDEvent *event)
{ // calculate distance between cluster and closest track

  Double_t trkPos[3] = {0,0,0};

  Int_t bestTrkMatchId = -1;
  Double_t distance = 9999; // init to a big number

  Double_t dist = 0;
  Double_t distX = 0, distY = 0, distZ = 0;
 
  for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
    AliESDtrack *track = event->GetTrack(iTrack);
    if (!track) {
      Printf("ERROR: Could not get track %d", iTrack);
      continue;
    }

    // check for approx. eta and phi range before we propagate..
    // TBD

    AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
    if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
      continue;
    }
    emctrack->GetXYZ(trkPos);
    if (emctrack) delete emctrack;

    distX = clsPos[0]-trkPos[0];
    distY = clsPos[1]-trkPos[1];
    distZ = clsPos[2]-trkPos[2];
    dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);

    if (dist < distance) {
      distance = dist;
      bestTrkMatchId = iTrack;
    }
  } // iTrack

  // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
  *trkMatchId = bestTrkMatchId;
  return distance;
}
Commit	Line	Data
cf6522d1	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Base class for ESD analysis
	4	// - reconstruction output
	5	// implementation file
	6	//
	7	//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
	8	//_________________________________________________________________________
2fbf38ac	9
	10	#include "AliAnalysisEtReconstructed.h"
	11	#include "AliAnalysisEtCuts.h"
	12	#include "AliESDtrack.h"
ba136eb4	13	#include "AliEMCALTrack.h"
2fbf38ac	14	#include "AliESDCaloCluster.h"
2fbf38ac	15	#include "TVector3.h"
ba136eb4	16	#include "TGeoGlobalMagField.h"
ba136eb4	17	#include "AliMagF.h"
2fbf38ac	18	#include "AliVEvent.h"
2fbf38ac	19	#include "AliESDEvent.h"
b5821c13	20	#include "AliESDtrackCuts.h"
2fbf38ac	21	#include "AliVParticle.h"
87efb15c	22	#include "TDatabasePDG.h"
87efb15c	23	#include "TList.h"
b5821c13	24	#include "AliESDpid.h"
2fbf38ac	25	#include <iostream>
	26	#include "TH2F.h"
	27
16abb579	28	using namespace std;
	29
	30	ClassImp(AliAnalysisEtReconstructed);
	31
	32
2fbf38ac	33	AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
2fbf38ac	34	AliAnalysisEt()
2fbf38ac	35	,fTrackDistanceCut(0)
87efb15c	36	,fPidCut(0)
2fbf38ac	37	,fClusterType(0)
87efb15c	38	,fHistChargedPionEnergyDeposit(0)
	39	,fHistProtonEnergyDeposit(0)
	40	,fHistAntiProtonEnergyDeposit(0)
	41	,fHistChargedKaonEnergyDeposit(0)
	42	,fHistMuonEnergyDeposit(0)
2fbf38ac	43	{
	44
	45	}
	46
cf6522d1	47	AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
2fbf38ac	48	{
cf6522d1	49	}
	50
	51	Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
	52	{ // analyse ESD event
2fbf38ac	53	ResetEventValues();
	54	AliESDEvent event = dynamic_cast<AliESDEvent>(ev);
	55
87efb15c	56	Double_t protonMass = fPdgDB->GetParticle("proton")->Mass();
87efb15c	57
b5821c13	58	//for PID
b5821c13	59	AliESDpid *pID = new AliESDpid();
6ad010c1	60	pID->MakePID(event);
6ad010c1	61	TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event);
b5821c13	62	Int_t nGoodTracks = list->GetEntries();
6ad010c1	63	// printf("nGoodTracks %d nCaloClusters %d\n", nGoodTracks, event->GetNumberOfCaloClusters());
6ad010c1	64
b5821c13	65	for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
	66	{
	67	AliESDtrack track = dynamic_cast<AliESDtrack> (list->At(iTrack));
2fbf38ac	68	if (!track)
	69	{
	70	Printf("ERROR: Could not get track %d", iTrack);
	71	continue;
	72	}
	73
	74	fMultiplicity++;
	75
b5821c13	76
	77	Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
	78	nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
	79	nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
	80	nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
	81	nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
16abb579	82	/*
b5821c13	83	bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
	84	bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
	85	bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
	86	bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
16abb579	87	*/
b5821c13	88
2fbf38ac	89	Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
	90	Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
	91
	92	Float_t massPart = 0;
	93
	94	const Double_t *pidWeights = track->PID();
87efb15c	95	Int_t maxpid = -1;
	96	Double_t maxpidweight = 0;
	97
2fbf38ac	98	if (pidWeights)
2fbf38ac	99	{
2fbf38ac	100	for (Int_t p =0; p < AliPID::kSPECIES; p++)
	101	{
	102	if (pidWeights[p] > maxpidweight)
	103	{
	104	maxpidweight = pidWeights[p];
	105	maxpid = p;
	106	}
	107	}
	108	if (maxpid == AliPID::kProton)
	109	{
87efb15c	110	//by definition of ET
87efb15c	111	massPart = -protonMass*track->Charge();
2fbf38ac	112	}
	113
	114	}
	115
	116	Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
6ad010c1	117	//printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout
2fbf38ac	118
4998becf	119	if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
2fbf38ac	120	{
	121	fTotChargedEt += et;
	122	fChargedMultiplicity++;
87efb15c	123	if (maxpid != -1)
	124	{
	125	if (maxpid == AliPID::kPion)
	126	{
	127	fProtonEt += et;
	128	}
	129	if (maxpid == AliPID::kKaon)
	130	{
	131	fChargedKaonEt += et;
	132	}
	133	if (maxpid == AliPID::kMuon)
	134	{
	135	fMuonEt += et;
	136	}
	137	if (maxpid == AliPID::kElectron)
	138	{
	139	fElectronEt += et;
	140	}
	141	}
2fbf38ac	142
	143	if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
	144	{
	145	fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
87efb15c	146	if (maxpid != -1)
	147	{
	148	if (maxpid == AliPID::kPion)
	149	{
	150	fProtonEtAcc += et;
	151	}
	152	if (maxpid == AliPID::kKaon)
	153	{
	154	fChargedKaonEtAcc += et;
	155	}
	156	if (maxpid == AliPID::kMuon)
	157	{
	158	fMuonEtAcc += et;
	159	}
	160	if (maxpid == AliPID::kElectron)
	161	{
	162	fElectronEtAcc += et;
	163	}
	164	}
	165
2fbf38ac	166	}
	167	}
	168
2fbf38ac	169	if (TrackHitsCalorimeter(track, event->GetMagneticField()))
2fbf38ac	170	{
13b0d3c1	171	Double_t phi = track->Phi();
13b0d3c1	172	Double_t pt = track->Pt();
cf6522d1	173	// printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
13b0d3c1	174	if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
13b0d3c1	175	else fHistPhivsPtNeg->Fill(phi, pt);
2fbf38ac	176	}
e6db391c	177	}
e6db391c	178	delete pID;
2fbf38ac	179
	180	for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
	181	{
	182	AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
	183	if (!cluster)
	184	{
	185	Printf("ERROR: Could not get cluster %d", iCluster);
	186	continue;
	187	}
	188
87efb15c	189	if (cluster->GetType() != fClusterType) continue;
ba136eb4	190	//if(cluster->GetTracksMatched() > 0)
6ad010c1	191	// printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout
87efb15c	192
87efb15c	193
2fbf38ac	194	if (cluster->E() < fClusterEnergyCut) continue;
	195	Float_t pos[3];
	196	TVector3 cp(pos);
	197	cluster->GetPosition(pos);
87efb15c	198
ba136eb4	199	Double_t distance = cluster->GetEmcCpvDistance();
	200	Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
	201	if ( cluster->IsEMCAL() ) {
	202	distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
	203	}
	204
	205	fHistTMDeltaR->Fill(distance);
	206	if (distance < fTrackDistanceCut)
2fbf38ac	207	{
ba136eb4	208	if (cluster->GetNTracksMatched() == 1 && trackMatchedIndex>0)
87efb15c	209	{
ba136eb4	210	AliVTrack *track = event->GetTrack(trackMatchedIndex);
87efb15c	211	const Double_t *pidWeights = track->PID();
	212
	213	Double_t maxpidweight = 0;
	214	Int_t maxpid = 0;
	215
	216	if (pidWeights)
	217	{
	218	for (Int_t p =0; p < AliPID::kSPECIES; p++)
	219	{
	220	if (pidWeights[p] > maxpidweight)
	221	{
	222	maxpidweight = pidWeights[p];
	223	maxpid = p;
	224	}
	225	}
	226	if(maxpidweight > fPidCut)
	227	{
	228	if(maxpid == AliPID::kProton)
	229	{
	230	if(track->Charge() == 1)
	231	{
	232	fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
	233	}
	234	else if(track->Charge() == -1)
	235	{
	236	fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
	237	}
	238	}
	239	else if(maxpid == AliPID::kPion)
	240	{
	241	fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
	242	}
	243	else if(maxpid == AliPID::kKaon)
	244	{
	245	fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
	246	}
	247	else if(maxpid == AliPID::kMuon)
	248	{
	249	fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
	250	}
	251	}
	252	}
	253	}
	254
2fbf38ac	255	continue;
ba136eb4	256	} // distance
2fbf38ac	257
4998becf	258	if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
2fbf38ac	259
	260	cluster->GetPosition(pos);
	261
	262	// TODO: replace with TVector3, too lazy now...
	263
	264	float dist = TMath::Sqrt(pos[0]pos[0] + pos[1]pos[1]);
	265
	266	float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
	267	// float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
	268	fTotNeutralEt += cluster->E() * TMath::Sin(theta);
13b0d3c1	269	fNeutralMultiplicity++;
2fbf38ac	270
2fbf38ac	271	fMultiplicity++;
2fbf38ac	272	}
	273
	274	fTotNeutralEtAcc = fTotNeutralEt;
	275	fTotEt = fTotChargedEt + fTotNeutralEt;
	276	fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
	277
2fbf38ac	278	// Fill the histograms...
	279	FillHistograms();
	280
	281	return 0;
	282	}
	283
	284	bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
cf6522d1	285	{ // check vertex
2fbf38ac	286
	287	Float_t bxy = 999.;
	288	Float_t bz = 999.;
	289	dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
	290
4998becf	291	bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
	292	(TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
	293	(TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
	294	(TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
	295	(TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
2fbf38ac	296
4998becf	297	return status;
2fbf38ac	298	}
	299
	300	void AliAnalysisEtReconstructed::Init()
cf6522d1	301	{ // Init
2fbf38ac	302	AliAnalysisEt::Init();
83d0f02c	303	fPidCut = fCuts->GetReconstructedPidCut();
ba136eb4	304	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
ba136eb4	305
2fbf38ac	306	}
	307
	308	bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
cf6522d1	309	{ // propagate track to detector radius
2fbf38ac	310
2fbf38ac	311	AliESDtrack esdTrack = dynamic_cast<AliESDtrack>(track);
cf6522d1	312	// Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac	313
	314	Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
	315
cf6522d1	316	// if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac	317	return prop &&
	318	TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
	319	esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
	320	esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
	321	}
	322
87efb15c	323	void AliAnalysisEtReconstructed::FillOutputList(TList* list)
ce546038	324	{ // add some extra histograms to the ones from base class
87efb15c	325	AliAnalysisEt::FillOutputList(list);
	326
	327	list->Add(fHistChargedPionEnergyDeposit);
	328	list->Add(fHistProtonEnergyDeposit);
	329	list->Add(fHistAntiProtonEnergyDeposit);
	330	list->Add(fHistChargedKaonEnergyDeposit);
	331	list->Add(fHistMuonEnergyDeposit);
	332	}
	333
	334	void AliAnalysisEtReconstructed::CreateHistograms()
ce546038	335	{ // add some extra histograms to the ones from base class
87efb15c	336	AliAnalysisEt::CreateHistograms();
87efb15c	337
0fa8c632	338	Int_t nbinsEt = 1000;
	339	Double_t minEt = 0;
	340	Double_t maxEt = 10;
	341
	342	// possibly change histogram limits
	343	if (fCuts) {
	344	nbinsEt = fCuts->GetHistNbinsParticleEt();
	345	minEt = fCuts->GetHistMinParticleEt();
	346	maxEt = fCuts->GetHistMaxParticleEt();
	347	}
	348
87efb15c	349	TString histname;
87efb15c	350	histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	351	fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	352	fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
	353	fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
	354
	355	histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	356	fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	357	fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
	358	fHistProtonEnergyDeposit->SetYTitle("Energy of track");
	359
	360	histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	361	fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	362	fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
	363	fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
	364
	365	histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	366	fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	367	fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
	368	fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
	369
	370	histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	371	fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	372	fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
	373	fHistMuonEnergyDeposit->SetYTitle("Energy of track");
	374
	375
	376	}
ba136eb4	377
ba136eb4	378	Double_t
ce546038	379	AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
ba136eb4	380	Int_t *trkMatchId,
ce546038	381	const AliESDEvent *event)
ba136eb4	382	{ // calculate distance between cluster and closest track
	383
	384	Double_t trkPos[3] = {0,0,0};
	385
	386	Int_t bestTrkMatchId = -1;
	387	Double_t distance = 9999; // init to a big number
	388
	389	Double_t dist = 0;
	390	Double_t distX = 0, distY = 0, distZ = 0;
	391
	392	for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
	393	AliESDtrack *track = event->GetTrack(iTrack);
	394	if (!track) {
	395	Printf("ERROR: Could not get track %d", iTrack);
	396	continue;
	397	}
	398
	399	// check for approx. eta and phi range before we propagate..
	400	// TBD
	401
	402	AliEMCALTrack emctrack = new AliEMCALTrack(track);
	403	if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
	404	continue;
	405	}
	406	emctrack->GetXYZ(trkPos);
	407	if (emctrack) delete emctrack;
	408
	409	distX = clsPos[0]-trkPos[0];
	410	distY = clsPos[1]-trkPos[1];
	411	distZ = clsPos[2]-trkPos[2];
	412	dist = TMath::Sqrt(distXdistX + distYdistY + distZ*distZ);
	413
	414	if (dist < distance) {
	415	distance = dist;
	416	bestTrkMatchId = iTrack;
	417	}
	418	} // iTrack
	419
6ad010c1	420	// printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
ba136eb4	421	*trkMatchId = bestTrkMatchId;
	422	return distance;
	423	}