[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 "AliESDCaloCluster.h"
#include "TVector3.h"
#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliVParticle.h"
#include <iostream>
#include "TH2F.h"

AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
        AliAnalysisEt()
        ,fNTpcClustersCut(EtReconstructedCuts::kNTpcClustersCut)
        ,fNItsClustersCut(EtReconstructedCuts::knItsClustersCut)
        ,fTrackDistanceCut(0)
        ,fClusterType(0)
{

}

AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed() 
{
}

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

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

        fMultiplicity++;

        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();
        if (pidWeights)
        {
            Int_t maxpid = -1;
            Float_t maxpidweight = 0;
            for (Int_t p =0; p < AliPID::kSPECIES; p++)
            {
                if (pidWeights[p] > maxpidweight)
                {
                    maxpidweight = pidWeights[p];
                    maxpid = p;
                }
            }
            if (maxpid == AliPID::kProton)
            {
                //     massPart = -0.938*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()) < fEtaCut && CheckGoodVertex(track) && nItsClusters > fNItsClustersCut && nTPCClusters > fNTpcClustersCut)
        {
	    fTotChargedEt +=  et;
            fChargedMultiplicity++;

            if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
            {
                fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
            }
        }

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

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

	// printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells()); // tmp/debug printout
        if (cluster->GetType() != fClusterType) continue;

        if (cluster->E() < fClusterEnergyCut) continue;
        Float_t pos[3];
        TVector3 cp(pos);
        cluster->GetPosition(pos);
        //if (pos[0] < -(32.0*2.2)) continue; //Ensure that modules 0 and 1 are not used
        // if(cp.Phi() < 260.*TMath::Pi()/180.) continue;
        fHistTMDeltaR->Fill(cluster->GetEmcCpvDistance());
        if (cluster->GetEmcCpvDistance() < fTrackDistanceCut)
        {
            continue;
            //AliVParticle *matchedTrack = event->GetTrack(cluster->GetTrackMatched());
// 	    if(CheckGoodVertex(matchedTrack))
// 	    {
// 	       totChargedEnergy +=  matchedTrack->E();;
// 	       totChargedEt += matchedTrack->E()*TMath::Sin(matchedTrack);
// 	    }
        }

        if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == EtCommonCuts::kSingleCell) 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;

    std::cout << fTotChargedEtAcc << std::endl;
    // 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);

    // printf("Rec CheckGoodVertex: TMath::Abs(track->Xv()) %f fVertexXCut %f TMath::Abs(track->Yv()) %f fVertexYCut %f TMath::Abs(track->Zv()) %f fVertexZCut %f TMath::Abs(bxy) %f fIPxyCut %f TMath::Abs(bz) %f fIPzCut %f\n", TMath::Abs(track->Xv()), fVertexXCut, TMath::Abs(track->Yv()), fVertexYCut, TMath::Abs(track->Zv()), fVertexZCut, TMath::Abs(bxy), fIPxyCut, TMath::Abs(bz), fIPzCut); // tmp/debug printout

    return TMath::Abs(track->Xv()) < fVertexXCut && TMath::Abs(track->Yv()) < fVertexYCut && TMath::Abs(track->Zv()) < fVertexZCut && TMath::Abs(bxy) < fIPxyCut && TMath::Abs(bz) < fIPzCut;

}

void AliAnalysisEtReconstructed::Init()
{ // Init
    AliAnalysisEt::Init();
}

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.;
}
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"
	13	#include "AliESDCaloCluster.h"
	14	#include "TVector3.h"
	15	#include "AliVEvent.h"
	16	#include "AliESDEvent.h"
	17	#include "AliVParticle.h"
	18	#include <iostream>
	19	#include "TH2F.h"
	20
	21	AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
	22	AliAnalysisEt()
	23	,fNTpcClustersCut(EtReconstructedCuts::kNTpcClustersCut)
	24	,fNItsClustersCut(EtReconstructedCuts::knItsClustersCut)
	25	,fTrackDistanceCut(0)
	26	,fClusterType(0)
	27	{
	28
	29	}
	30
cf6522d1	31	AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
2fbf38ac	32	{
cf6522d1	33	}
	34
	35	Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
	36	{ // analyse ESD event
2fbf38ac	37	ResetEventValues();
	38	AliESDEvent event = dynamic_cast<AliESDEvent>(ev);
	39
	40	for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++)
	41	{
	42	AliVParticle *track = event->GetTrack(iTrack);
	43	if (!track)
	44	{
	45	Printf("ERROR: Could not get track %d", iTrack);
	46	continue;
	47	}
	48
	49	fMultiplicity++;
	50
	51	Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
	52	Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
	53
	54	Float_t massPart = 0;
	55
	56	const Double_t *pidWeights = track->PID();
	57	if (pidWeights)
	58	{
	59	Int_t maxpid = -1;
	60	Float_t maxpidweight = 0;
	61	for (Int_t p =0; p < AliPID::kSPECIES; p++)
	62	{
	63	if (pidWeights[p] > maxpidweight)
	64	{
	65	maxpidweight = pidWeights[p];
	66	maxpid = p;
	67	}
	68	}
	69	if (maxpid == AliPID::kProton)
	70	{
	71	// massPart = -0.938*track->Charge();
	72	}
	73
	74	}
	75
	76	Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
cf6522d1	77	// 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	78
	79	if (TMath::Abs(track->Eta()) < fEtaCut && CheckGoodVertex(track) && nItsClusters > fNItsClustersCut && nTPCClusters > fNTpcClustersCut)
	80	{
	81	fTotChargedEt += et;
	82	fChargedMultiplicity++;
	83
	84	if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
	85	{
	86	fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
	87	}
	88	}
	89
2fbf38ac	90	if (TrackHitsCalorimeter(track, event->GetMagneticField()))
2fbf38ac	91	{
13b0d3c1	92	Double_t phi = track->Phi();
13b0d3c1	93	Double_t pt = track->Pt();
cf6522d1	94	// printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
13b0d3c1	95	if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
13b0d3c1	96	else fHistPhivsPtNeg->Fill(phi, pt);
2fbf38ac	97	}
	98	}
	99
	100	for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
	101	{
	102	AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
	103	if (!cluster)
	104	{
	105	Printf("ERROR: Could not get cluster %d", iCluster);
	106	continue;
	107	}
	108
cf6522d1	109	// printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells()); // tmp/debug printout
2fbf38ac	110	if (cluster->GetType() != fClusterType) continue;
	111
	112	if (cluster->E() < fClusterEnergyCut) continue;
	113	Float_t pos[3];
	114	TVector3 cp(pos);
	115	cluster->GetPosition(pos);
	116	//if (pos[0] < -(32.0*2.2)) continue; //Ensure that modules 0 and 1 are not used
	117	// if(cp.Phi() < 260.*TMath::Pi()/180.) continue;
	118	fHistTMDeltaR->Fill(cluster->GetEmcCpvDistance());
	119	if (cluster->GetEmcCpvDistance() < fTrackDistanceCut)
	120	{
	121	continue;
	122	//AliVParticle *matchedTrack = event->GetTrack(cluster->GetTrackMatched());
	123	// if(CheckGoodVertex(matchedTrack))
	124	// {
	125	// totChargedEnergy += matchedTrack->E();;
	126	// totChargedEt += matchedTrack->E()*TMath::Sin(matchedTrack);
	127	// }
	128	}
	129
	130	if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == EtCommonCuts::kSingleCell) continue;
	131
	132	cluster->GetPosition(pos);
	133
	134	// TODO: replace with TVector3, too lazy now...
	135
	136	float dist = TMath::Sqrt(pos[0]pos[0] + pos[1]pos[1]);
	137
	138	float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
	139	// float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
	140	fTotNeutralEt += cluster->E() * TMath::Sin(theta);
13b0d3c1	141	fNeutralMultiplicity++;
2fbf38ac	142
2fbf38ac	143	fMultiplicity++;
2fbf38ac	144	}
	145
	146	fTotNeutralEtAcc = fTotNeutralEt;
13b0d3c1	147
2fbf38ac	148	fTotEt = fTotChargedEt + fTotNeutralEt;
	149	fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
	150
	151	std::cout << fTotChargedEtAcc << std::endl;
	152	// Fill the histograms...
	153	FillHistograms();
	154
	155	return 0;
	156	}
	157
	158	bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
cf6522d1	159	{ // check vertex
2fbf38ac	160
	161	Float_t bxy = 999.;
	162	Float_t bz = 999.;
	163	dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
	164
cf6522d1	165	// printf("Rec CheckGoodVertex: TMath::Abs(track->Xv()) %f fVertexXCut %f TMath::Abs(track->Yv()) %f fVertexYCut %f TMath::Abs(track->Zv()) %f fVertexZCut %f TMath::Abs(bxy) %f fIPxyCut %f TMath::Abs(bz) %f fIPzCut %f\n", TMath::Abs(track->Xv()), fVertexXCut, TMath::Abs(track->Yv()), fVertexYCut, TMath::Abs(track->Zv()), fVertexZCut, TMath::Abs(bxy), fIPxyCut, TMath::Abs(bz), fIPzCut); // tmp/debug printout
	166
	167	return TMath::Abs(track->Xv()) < fVertexXCut && TMath::Abs(track->Yv()) < fVertexYCut && TMath::Abs(track->Zv()) < fVertexZCut && TMath::Abs(bxy) < fIPxyCut && TMath::Abs(bz) < fIPzCut;
2fbf38ac	168
	169	}
	170
	171	void AliAnalysisEtReconstructed::Init()
cf6522d1	172	{ // Init
2fbf38ac	173	AliAnalysisEt::Init();
2fbf38ac	174	}
	175
	176	bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
cf6522d1	177	{ // propagate track to detector radius
2fbf38ac	178
2fbf38ac	179	AliESDtrack esdTrack = dynamic_cast<AliESDtrack>(track);
cf6522d1	180	// Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac	181
	182	Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
	183
cf6522d1	184	// if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac	185	return prop &&
	186	TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
	187	esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
	188	esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
	189	}
	190