[u/mrichter/AliRoot.git] / PWGLF / 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"
#include "TH2I.h"
#include "TH1I.h"
#include "TFile.h"
#include "AliAnalysisHadEtCorrections.h"
#include "AliAnalysisEtSelector.h"
#include "AliLog.h"
#include "AliCentrality.h"
#include "AliPHOSGeoUtils.h"
#include "AliPHOSGeometry.h"


using namespace std;

ClassImp(AliAnalysisEtReconstructed);


AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
        AliAnalysisEt()
        ,fCorrections(0)
        ,fPidCut(0)
        ,fHistChargedPionEnergyDeposit(0)
        ,fHistProtonEnergyDeposit(0)
        ,fHistAntiProtonEnergyDeposit(0)
        ,fHistChargedKaonEnergyDeposit(0)
        ,fHistMuonEnergyDeposit(0)
        ,fHistRemovedEnergy(0)
        ,fGeomCorrection(1.0)
        ,fEMinCorrection(1.0/0.89)
	,fRecEffCorrection(1.0)
	,fClusterPosition(0)
	,fHistChargedEnergyRemoved(0)
	,fHistNeutralEnergyRemoved(0)
	,fHistGammaEnergyAdded(0)
{

}

AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
{//destructor
    delete fCorrections;
    delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
    delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
    delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
    delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
    delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */

    delete fHistRemovedEnergy; // removed energy
    delete fClusterPosition;
    delete fHistChargedEnergyRemoved;
    delete fHistNeutralEnergyRemoved;
    delete fHistGammaEnergyAdded;

}

Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
{

    //AliAnalysisEt::AnalyseEvent(ev);
    // analyse ESD event
    ResetEventValues();
    if (!ev) {
        AliFatal("ERROR: Event does not exist");
        return 0;
    }

    AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
    if (!event) {
        AliFatal("ERROR: ESD Event does not exist");
        return 0;
    }
    fSelector->SetEvent(event);
    
    Int_t cent = -1;
    if (fCentrality)
    {
        cent = fCentrality->GetCentralityClass10("V0M");
        fCentClass = fCentrality->GetCentralityClass10("V0M");
    }

    for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
    {
        AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
        if (!cluster)
        {
            AliError(Form("ERROR: Could not get cluster %d", iCluster));
            continue;
        }
        int x = 0;
	fCutFlow->Fill(x++);
	if(cluster->IsEMCAL()) continue;
	fCutFlow->Fill(x++);
	if(!fSelector->CutMinEnergy(*cluster)) continue;
	fCutFlow->Fill(x++);
        if (!fSelector->CutDistanceToBadChannel(*cluster)) continue;
	fCutFlow->Fill(x++);

        Float_t pos[3];

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

        Bool_t matched = false;

	
	matched = !fSelector->CutTrackMatching(*cluster);

        if (matched)
        {
	  
            if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
            {
                AliVTrack *track = event->GetTrack(trackMatchedIndex);
                if (!track) {
                    AliError("Error: track does not exist");
                }
                else {
                    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 (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
                        {
                            fEnergyDeposited = cluster->E();
                            fMomentumTPC = track->P();
                            fCharge = track->Charge();
                            fParticlePid = maxpid;
                            fPidProb = maxpidweight;
                            AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
                            if (!esdTrack) {
                                AliError("Error: track does not exist");
                            }
                            else {
                                if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
                                fDepositTree->Fill();
                            }
                        }

                        if (maxpidweight > fPidCut)
                        {
                            //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);

                            //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;

                            //Float_t et = cluster->E() * TMath::Sin(theta);
                            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
        else
        {
	  fCutFlow->Fill(x++);
	  //std::cout << x++ << std::endl;

            //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
            //if (cluster->E() < fClusterEnergyCut) continue;
            cluster->GetPosition(pos);
	    
	    TVector3 p2(pos);
	    
	    fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());

	    fTotNeutralEt += CalculateTransverseEnergy(cluster);
            fNeutralMultiplicity++;
        }
        fMultiplicity++;
    }
    
    fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
    fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
    fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
    fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
    
    fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
    fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);

    Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
    fHistRemovedEnergy->Fill(removedEnergy);
    
    fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
    fTotEt = fTotChargedEt + fTotNeutralEt;
// Fill the histograms...0
    FillHistograms();
   // std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
    return 0;
}

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

    Float_t bxy = 999.;
    Float_t bz = 999.;
    if (!track) {
        AliError("ERROR: no track");
        return kFALSE;
    }
    AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
    if (!esdTrack) {
        AliError("ERROR: no track");
        return kFALSE;
    }
    esdTrack->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));
    if (!fCorrections) {
        cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
    }
}

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

    if (!track) {
        cout<<"Warning: track empty"<<endl;
        return kFALSE;
    }
    AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
    if (!esdTrack) {
        AliError("ERROR: no ESD track");
        return kFALSE;
    }
    // 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 && fSelector->CutGeometricalAcceptance(*esdTrack);
}

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);

    list->Add(fHistRemovedEnergy);
    list->Add(fClusterPosition);
    
    list->Add(fHistChargedEnergyRemoved);
    list->Add(fHistNeutralEnergyRemoved);
    list->Add(fHistGammaEnergyAdded);
}

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");

    histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
    fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
    //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
    //fHistMuonEnergyDeposit->SetYTitle("Energy of track");

    histname = "fClusterPosition" + fHistogramNameSuffix;
    fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13);
    fClusterPosition->SetXTitle("Energy deposited in calorimeter");
    fClusterPosition->SetYTitle("Energy of track");


    histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
    fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);

    histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
    fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);

    histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
    fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);


}

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) {
            AliError(Form("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);
        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>
2fbf38ac	26	#include "TH2F.h"
ef647350	27	#include "TH2I.h"
	28	#include "TH1I.h"
	29	#include "TFile.h"
964c8159	30	#include "AliAnalysisHadEtCorrections.h"
ef647350	31	#include "AliAnalysisEtSelector.h"
0f6416f3	32	#include "AliLog.h"
e9da35da	33	#include "AliCentrality.h"
ef647350	34	#include "AliPHOSGeoUtils.h"
ef647350	35	#include "AliPHOSGeometry.h"
0f6416f3	36
2fbf38ac	37
16abb579	38	using namespace std;
	39
	40	ClassImp(AliAnalysisEtReconstructed);
	41
	42
2fbf38ac	43	AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
f61cec2f	44	AliAnalysisEt()
	45	,fCorrections(0)
	46	,fPidCut(0)
	47	,fHistChargedPionEnergyDeposit(0)
	48	,fHistProtonEnergyDeposit(0)
	49	,fHistAntiProtonEnergyDeposit(0)
	50	,fHistChargedKaonEnergyDeposit(0)
	51	,fHistMuonEnergyDeposit(0)
	52	,fHistRemovedEnergy(0)
	53	,fGeomCorrection(1.0)
	54	,fEMinCorrection(1.0/0.89)
	55	,fRecEffCorrection(1.0)
	56	,fClusterPosition(0)
	57	,fHistChargedEnergyRemoved(0)
	58	,fHistNeutralEnergyRemoved(0)
	59	,fHistGammaEnergyAdded(0)
2fbf38ac	60	{
	61
	62	}
	63
e9da35da	64	AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
f61cec2f	65	{//destructor
e9da35da	66	delete fCorrections;
ef647350	67	delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
	68	delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
	69	delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
	70	delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
d0c22dcc	71	delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
	72
	73	delete fHistRemovedEnergy; // removed energy
311c6540	74	delete fClusterPosition;
	75	delete fHistChargedEnergyRemoved;
	76	delete fHistNeutralEnergyRemoved;
	77	delete fHistGammaEnergyAdded;
d0c22dcc	78
cf6522d1	79	}
	80
	81	Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
e9da35da	82	{
f61cec2f	83
f61cec2f	84	//AliAnalysisEt::AnalyseEvent(ev);
e9da35da	85	// analyse ESD event
2fbf38ac	86	ResetEventValues();
e9da35da	87	if (!ev) {
	88	AliFatal("ERROR: Event does not exist");
	89	return 0;
	90	}
	91
2fbf38ac	92	AliESDEvent event = dynamic_cast<AliESDEvent>(ev);
e9da35da	93	if (!event) {
	94	AliFatal("ERROR: ESD Event does not exist");
	95	return 0;
ec956c46	96	}
ef647350	97	fSelector->SetEvent(event);
f61cec2f	98
743ce29f	99	Int_t cent = -1;
	100	if (fCentrality)
	101	{
	102	cent = fCentrality->GetCentralityClass10("V0M");
ef647350	103	fCentClass = fCentrality->GetCentralityClass10("V0M");
743ce29f	104	}
743ce29f	105
2fbf38ac	106	for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
	107	{
	108	AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
	109	if (!cluster)
	110	{
e9da35da	111	AliError(Form("ERROR: Could not get cluster %d", iCluster));
2fbf38ac	112	continue;
2fbf38ac	113	}
ef647350	114	int x = 0;
f61cec2f	115	fCutFlow->Fill(x++);
	116	if(cluster->IsEMCAL()) continue;
	117	fCutFlow->Fill(x++);
	118	if(!fSelector->CutMinEnergy(*cluster)) continue;
	119	fCutFlow->Fill(x++);
ef647350	120	if (!fSelector->CutDistanceToBadChannel(*cluster)) continue;
f61cec2f	121	fCutFlow->Fill(x++);
e9da35da	122
2fbf38ac	123	Float_t pos[3];
e9da35da	124
2fbf38ac	125	cluster->GetPosition(pos);
e9da35da	126	TVector3 cp(pos);
f61cec2f	127
f61cec2f	128	Double_t distance = cluster->GetEmcCpvDistance();
e9da35da	129	Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
	130	if ( cluster->IsEMCAL() ) {
	131	distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
	132	}
87efb15c	133
e9da35da	134	Bool_t matched = false;
e9da35da	135
f61cec2f	136
f61cec2f	137	matched = !fSelector->CutTrackMatching(*cluster);
e9da35da	138
	139	if (matched)
	140	{
f61cec2f	141
e9da35da	142	if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
	143	{
	144	AliVTrack *track = event->GetTrack(trackMatchedIndex);
	145	if (!track) {
	146	AliError("Error: track does not exist");
	147	}
	148	else {
	149	const Double_t *pidWeights = track->PID();
	150
	151	Double_t maxpidweight = 0;
	152	Int_t maxpid = 0;
	153
	154	if (pidWeights)
	155	{
	156	for (Int_t p =0; p < AliPID::kSPECIES; p++)
	157	{
	158	if (pidWeights[p] > maxpidweight)
	159	{
	160	maxpidweight = pidWeights[p];
	161	maxpid = p;
	162	}
	163	}
f61cec2f	164	if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
e9da35da	165	{
e9da35da	166	fEnergyDeposited = cluster->E();
f61cec2f	167	fMomentumTPC = track->P();
e9da35da	168	fCharge = track->Charge();
	169	fParticlePid = maxpid;
	170	fPidProb = maxpidweight;
	171	AliESDtrack esdTrack = dynamic_cast<AliESDtrack>(track);
	172	if (!esdTrack) {
	173	AliError("Error: track does not exist");
	174	}
	175	else {
	176	if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
f61cec2f	177	fDepositTree->Fill();
e9da35da	178	}
	179	}
	180
	181	if (maxpidweight > fPidCut)
	182	{
f61cec2f	183	//Float_t dist = TMath::Sqrt(pos[0]pos[0] + pos[1]pos[1]);
e9da35da	184
f61cec2f	185	//Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
e9da35da	186
f61cec2f	187	//Float_t et = cluster->E() * TMath::Sin(theta);
e9da35da	188	if (maxpid == AliPID::kProton)
	189	{
	190
	191	if (track->Charge() == 1)
	192	{
e9da35da	193	fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
	194	}
	195	else if (track->Charge() == -1)
	196	{
e9da35da	197	fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
	198	}
	199	}
	200	else if (maxpid == AliPID::kPion)
	201	{
e9da35da	202	fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
	203	}
	204	else if (maxpid == AliPID::kKaon)
	205	{
e9da35da	206	fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
	207	}
	208	else if (maxpid == AliPID::kMuon)
	209	{
	210	fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
	211	}
	212	}
	213	}
	214	}
	215	}
	216	//continue;
ba136eb4	217	} // distance
e9da35da	218	else
e9da35da	219	{
f61cec2f	220	fCutFlow->Fill(x++);
f61cec2f	221	//std::cout << x++ << std::endl;
e9da35da	222
ef647350	223	//if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
ef647350	224	//if (cluster->E() < fClusterEnergyCut) continue;
e9da35da	225	cluster->GetPosition(pos);
f61cec2f	226
	227	TVector3 p2(pos);
	228
	229	fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
2fbf38ac	230
f61cec2f	231	fTotNeutralEt += CalculateTransverseEnergy(cluster);
e9da35da	232	fNeutralMultiplicity++;
e9da35da	233	}
ef647350	234	fMultiplicity++;
ef647350	235	}
f61cec2f	236
ef647350	237	fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
	238	fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
	239	fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
	240	fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
f61cec2f	241
ef647350	242	fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
ef647350	243	fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
2fbf38ac	244
f61cec2f	245	Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
e9da35da	246	fHistRemovedEnergy->Fill(removedEnergy);
f61cec2f	247
ef647350	248	fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
2fbf38ac	249	fTotEt = fTotChargedEt + fTotNeutralEt;
f61cec2f	250	// Fill the histograms...0
2fbf38ac	251	FillHistograms();
f61cec2f	252	// std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
2fbf38ac	253	return 0;
	254	}
	255
	256	bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
f61cec2f	257	{ // check vertex
2fbf38ac	258
	259	Float_t bxy = 999.;
	260	Float_t bz = 999.;
e9da35da	261	if (!track) {
	262	AliError("ERROR: no track");
	263	return kFALSE;
ec956c46	264	}
1b8c3d66	265	AliESDtrack esdTrack = dynamic_cast<AliESDtrack>(track);
e9da35da	266	if (!esdTrack) {
	267	AliError("ERROR: no track");
	268	return kFALSE;
1b8c3d66	269	}
1b8c3d66	270	esdTrack->GetImpactParametersTPC(bxy,bz);
ec956c46	271
2fbf38ac	272
e9da35da	273	bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
	274	(TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
	275	(TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
	276	(TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
	277	(TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
2fbf38ac	278
4998becf	279	return status;
2fbf38ac	280	}
	281
	282	void AliAnalysisEtReconstructed::Init()
f61cec2f	283	{ // Init
2fbf38ac	284	AliAnalysisEt::Init();
83d0f02c	285	fPidCut = fCuts->GetReconstructedPidCut();
ba136eb4	286	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
e9da35da	287	if (!fCorrections) {
e9da35da	288	cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl;
0f97be4c	289	}
2fbf38ac	290	}
	291
	292	bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
f61cec2f	293	{ // propagate track to detector radius
2fbf38ac	294
e9da35da	295	if (!track) {
	296	cout<<"Warning: track empty"<<endl;
	297	return kFALSE;
	298	}
	299	AliESDtrack esdTrack= dynamic_cast<AliESDtrack>(track);
	300	if (!esdTrack) {
	301	AliError("ERROR: no ESD track");
	302	return kFALSE;
	303	}
	304	// Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac	305
	306	Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
	307
cf6522d1	308	// if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
f61cec2f	309	return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
2fbf38ac	310	}
2fbf38ac	311
87efb15c	312	void AliAnalysisEtReconstructed::FillOutputList(TList* list)
f61cec2f	313	{ // add some extra histograms to the ones from base class
87efb15c	314	AliAnalysisEt::FillOutputList(list);
	315
	316	list->Add(fHistChargedPionEnergyDeposit);
	317	list->Add(fHistProtonEnergyDeposit);
	318	list->Add(fHistAntiProtonEnergyDeposit);
	319	list->Add(fHistChargedKaonEnergyDeposit);
	320	list->Add(fHistMuonEnergyDeposit);
e9da35da	321
e9da35da	322	list->Add(fHistRemovedEnergy);
ef647350	323	list->Add(fClusterPosition);
f61cec2f	324
ef647350	325	list->Add(fHistChargedEnergyRemoved);
	326	list->Add(fHistNeutralEnergyRemoved);
	327	list->Add(fHistGammaEnergyAdded);
87efb15c	328	}
	329
	330	void AliAnalysisEtReconstructed::CreateHistograms()
f61cec2f	331	{ // add some extra histograms to the ones from base class
87efb15c	332	AliAnalysisEt::CreateHistograms();
87efb15c	333
0fa8c632	334	Int_t nbinsEt = 1000;
	335	Double_t minEt = 0;
	336	Double_t maxEt = 10;
	337
	338	// possibly change histogram limits
	339	if (fCuts) {
e9da35da	340	nbinsEt = fCuts->GetHistNbinsParticleEt();
	341	minEt = fCuts->GetHistMinParticleEt();
	342	maxEt = fCuts->GetHistMaxParticleEt();
0fa8c632	343	}
0fa8c632	344
87efb15c	345	TString histname;
87efb15c	346	histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	347	fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	348	fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
87efb15c	349	fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
e9da35da	350
87efb15c	351	histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	352	fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	353	fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
87efb15c	354	fHistProtonEnergyDeposit->SetYTitle("Energy of track");
e9da35da	355
87efb15c	356	histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	357	fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	358	fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
87efb15c	359	fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
e9da35da	360
87efb15c	361	histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	362	fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	363	fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
87efb15c	364	fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
e9da35da	365
87efb15c	366	histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632	367	fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c	368	fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
87efb15c	369	fHistMuonEnergyDeposit->SetYTitle("Energy of track");
e9da35da	370
	371	histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
	372	fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
	373	//fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
	374	//fHistMuonEnergyDeposit->SetYTitle("Energy of track");
	375
ef647350	376	histname = "fClusterPosition" + fHistogramNameSuffix;
	377	fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13);
	378	fClusterPosition->SetXTitle("Energy deposited in calorimeter");
	379	fClusterPosition->SetYTitle("Energy of track");
	380
	381
	382	histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
	383	fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
	384
	385	histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
	386	fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
	387
	388	histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
	389	fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
e9da35da	390
87efb15c	391
87efb15c	392	}
ba136eb4	393
e9da35da	394	Double_t
ce546038	395	AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
e9da35da	396	Int_t *trkMatchId,
e9da35da	397	const AliESDEvent *event)
f61cec2f	398	{ // calculate distance between cluster and closest track
ba136eb4	399
e9da35da	400	Double_t trkPos[3] = {0,0,0};
ba136eb4	401
e9da35da	402	Int_t bestTrkMatchId = -1;
e9da35da	403	Double_t distance = 9999; // init to a big number
ba136eb4	404
e9da35da	405	Double_t dist = 0;
e9da35da	406	Double_t distX = 0, distY = 0, distZ = 0;
ba136eb4	407
e9da35da	408	for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
	409	AliESDtrack *track = event->GetTrack(iTrack);
	410	if (!track) {
	411	AliError(Form("ERROR: Could not get track %d", iTrack));
	412	continue;
	413	}
ba136eb4	414
e9da35da	415	// check for approx. eta and phi range before we propagate..
e9da35da	416	// TBD
ba136eb4	417
e9da35da	418	AliEMCALTrack emctrack = new AliEMCALTrack(track);
	419	if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
	420	continue;
	421	}
	422	emctrack->GetXYZ(trkPos);
	423	delete emctrack;
ba136eb4	424
e9da35da	425	distX = clsPos[0]-trkPos[0];
	426	distY = clsPos[1]-trkPos[1];
	427	distZ = clsPos[2]-trkPos[2];
	428	dist = TMath::Sqrt(distXdistX + distYdistY + distZ*distZ);
ba136eb4	429
e9da35da	430	if (dist < distance) {
	431	distance = dist;
	432	bestTrkMatchId = iTrack;
	433	}
	434	} // iTrack
	435
	436	// printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
	437	*trkMatchId = bestTrkMatchId;
	438	return distance;
ba136eb4	439	}