- adding histograms for energy deposited by identified charged particles

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

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

#include "AliAnalysisEtMonteCarlo.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
#include "AliStack.h"
#include "AliMCEvent.h"
#include "TH2F.h"
#include "TParticle.h"

Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
{ // analyse MC event
     ResetEventValues();
     
    // Get us an mc event
    AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);

    // Let's play with the stack!
    AliStack *stack = event->Stack();

    Int_t nPrim = stack->GetNtrack();

    Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle

    for (Int_t iPart = 0; iPart < nPrim; iPart++)
    {

        TParticle *part = stack->Particle(iPart);

        if (!part)
        {
            Printf("ERROR: Could not get particle %d", iPart);
            continue;
        }

        TParticlePDG *pc = part->GetPDG(0);

        // Check if it is a primary particle
        if (!stack->IsPhysicalPrimary(iPart)) continue;

	// printf("MC: iPart %03d eta %4.3f phi %4.3f code %d charge %g \n", iPart, part->Eta(), part->Phi(), pc->PdgCode(), pc->Charge()); // tmp/debug printout

        // Check for reasonable (for now neutral and singly charged) charge on the particle
        //TODO:Maybe not only singly charged?
        if (TMath::Abs(pc->Charge()) != fCuts->GetMonteCarloSingleChargedParticle() && pc->Charge() != fCuts->GetMonteCarloNeutralParticle()) continue;

        fMultiplicity++;

        if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
        {

	  TParticlePDG *pdgCode =  part->GetPDG(0);
            if (
                TMath::Abs(pdgCode->PdgCode()) == fProtonCode ||
                TMath::Abs(pdgCode->PdgCode()) == fNeutronCode ||
                TMath::Abs(pdgCode->PdgCode()) == fLambdaCode ||
                TMath::Abs(pdgCode->PdgCode()) == fXiCode ||
                TMath::Abs(pdgCode->PdgCode()) == fXi0Code ||
                TMath::Abs(pdgCode->PdgCode()) == fOmegaCode
	       )
            {
                particleMassPart = -TMath::Sign(pdgCode->PdgCode(), pdgCode->PdgCode())*pdgCode->Mass();
            }
            
            if (pdgCode->Charge() == fCuts->GetMonteCarloNeutralParticle() )
            {
	       fNeutralMultiplicity++;
                fTotNeutralEt += part->Energy()*TMath::Sin(part->Theta());

                if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
                {
                    fTotNeutralEtAcc += part->Energy()*TMath::Sin(part->Theta());
                    fTotEtAcc += part->Energy()*TMath::Sin(part->Theta());
                }
            }
            else if (pdgCode->Charge() != fCuts->GetMonteCarloNeutralParticle() )
            {
	       fChargedMultiplicity++;
                fTotChargedEt += part->Energy()*TMath::Sin(part->Theta());
                if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
                {
                    fTotChargedEtAcc += part->Energy()*TMath::Sin(part->Theta());
                    fTotEtAcc += part->Energy()*TMath::Sin(part->Theta());
                }

		//	  if (TrackHitsCalorimeter(part, event->GetMagneticField()))
		if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
		  {
		    if (pdgCode->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
		    else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
		  }
	    }
	}
    }
    
    fTotNeutralEtAcc = fTotNeutralEt;
    fTotEt = fTotChargedEt + fTotNeutralEt;
    fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
    
    FillHistograms();

    return 0;    
}

void AliAnalysisEtMonteCarlo::Init()
{
    AliAnalysisEt::Init();
}

bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
{
  //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
   AliESDtrack *esdTrack = new AliESDtrack(part);
   // Printf("MC 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());

    bool status = prop && 
      TMath::Abs(esdTrack->Eta()) < fEtaCutAcc && 
      esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. && 
      esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
    delete esdTrack;

    return status;
}