[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.687)
	,fRecEffCorrection(1.0)
	,fClusterPosition(0)
	,fClusterEnergy(0)
	,fClusterEt(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 fClusterEnergy;
    delete fClusterEt;
    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;
    }
    if(!fSelector){
        AliFatal("ERROR: fSelector does not exist");
        return 0;
    }
    fSelector->SetEvent(event);
    
    Int_t cent = -1;
    if (fCentrality && cent)
    {
        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(!fSelector->IsDetectorCluster(*cluster)) continue;
	fCutFlow->Fill(x++);
	if(!fSelector->PassMinEnergyCut(*cluster)) continue;
	fCutFlow->Fill(x++);
        if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
	fCutFlow->Fill(x++);

        Float_t pos[3];

        cluster->GetPosition(pos);
        TVector3 cp(pos);
	
        Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();

        Bool_t matched = false;

	
	matched = !fSelector->PassTrackMatchingCut(*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());
	    fClusterEnergy->Fill(cluster->E());
	    fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E());

	    fTotNeutralEt += CorrectForReconstructionEfficiency(*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(fClusterEnergy);
    list->Add(fClusterEt);
    
    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 = "fClusterEnergy" + fHistogramNameSuffix;
    fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
    fClusterEnergy->SetXTitle("Number of clusters");
    fClusterEnergy->SetYTitle("Energy of cluster");

    histname = "fClusterEt" + fHistogramNameSuffix;
    fClusterEt = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
    fClusterEt->SetXTitle("Number of clusters");
    fClusterEt->SetYTitle("E_{T} of cluster");

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


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