<commit_before>//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Base class for MC analysis
//  - MC output
//  implementation file
//
//*-- Authors: <NAME> (Bergen), <NAME> (ORNL)
//_________________________________________________________________________

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

using namespace std;

ClassImp(AliAnalysisEtMonteCarlo);


// ctor
AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo() :
  AliAnalysisEt()
  ,fImpactParameter(0)
  ,fNcoll(0)
  ,fNpart(0)
{
}

// dtor
AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo() 
{
}

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

    // Hijing header
    AliGenEventHeader* genHeader = event->GenEventHeader();
    AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
    if (hijingGenHeader) {
      fImpactParameter = hijingGenHeader->ImpactParameter();
      fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
      fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants(); 
      /*
      printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
	     hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
      printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
	     hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants()); 
      printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
	     hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
      printf("NN %d NNw %d NwN %d, NwNw %d\n",
	     hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
      */
    }

    /* // placeholder if we want to get some Pythia info later
    AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
    if (pythiaGenHeader) { // not Hijing; try with Pythia      
      printf("Pythia: ProcessType %d  GetPtHard %g \n",
	     pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
    }
    */

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

    Int_t nPrim = stack->GetNtrack();

    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 *pdg = part->GetPDG(0);

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

        // 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(), pdg->PdgCode(), pdg->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(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;

        fMultiplicity++;

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

            if (
                TMath::Abs(pdg->PdgCode()) == fProtonCode ||
                TMath::Abs(pdg->PdgCode()) == fNeutronCode ||
                TMath::Abs(pdg->PdgCode()) == fLambdaCode ||
                TMath::Abs(pdg->PdgCode()) == fXiCode ||
                TMath::Abs(pdg->PdgCode()) == fXi0Code ||
                TMath::Abs(pdg->PdgCode()) == fOmegaCode
	       )
            {
	      if (pdg->PdgCode() > 0) { particleMassPart = - pdg->Mass();}
	      if (pdg->PdgCode() < 0) { particleMassPart = pdg->Mass();}
	    }
	    Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
	    	      
	    if (pdg->PdgCode() == fProtonCode || pdg->PdgCode() == fAntiProtonCode)
	      {
		fProtonEt += et;
	      }
	    if (pdg->PdgCode() == fPiPlusCode || pdg->PdgCode() == fPiMinusCode)
	      {
		fPionEt += et;
	      }
	    if (pdg->PdgCode() == fKPlusCode || pdg->PdgCode() == fKMinusCode)
	      {
		fChargedKaonEt += et;
	      }
	    if (pdg->PdgCode() == fMuPlusCode || pdg->PdgCode() == fMuMinusCode)
	      {
		fMuonEt += et;
	      }
	    if (pdg->PdgCode() == fEPlusCode || pdg->PdgCode() == fEMinusCode)<commit_msg>correcting baryon mass subtraction for visible energy in MC<commit_after>//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Base class for MC analysis
//  - MC output
//  implementation file
//
//*-- Authors: <NAME> (Bergen), <NAME> (ORNL)
//_________________________________________________________________________

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

using namespace std;

ClassImp(AliAnalysisEtMonteCarlo);


// ctor
AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo() :
  AliAnalysisEt()
  ,fImpactParameter(0)
  ,fNcoll(0)
  ,fNpart(0)
{
}

// dtor
AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo() 
{
}

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

    Double_t protonMass = fPdgDB->GetParticle("proton")->Mass(); // should maybe use average of proton and neutron, but they are pretty close

    // Hijing header
    AliGenEventHeader* genHeader = event->GenEventHeader();
    AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
    if (hijingGenHeader) {
      fImpactParameter = hijingGenHeader->ImpactParameter();
      fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
      fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants(); 
      /*
      printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
	     hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
      printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
	     hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants()); 
      printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
	     hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
      printf("NN %d NNw %d NwN %d, NwNw %d\n",
	     hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
      */
    }

    /* // placeholder if we want to get some Pythia info later
    AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
    if (pythiaGenHeader) { // not Hijing; try with Pythia      
      printf("Pythia: ProcessType %d  GetPtHard %g \n",
	     pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
    }
    */

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

    Int_t nPrim = stack->GetNtrack();

    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 *pdg = part->GetPDG(0);

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

        // 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(), pdg->PdgCode(), pdg->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(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;

        fMultiplicity++;

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

            if (
                TMath::Abs(pdg->PdgCode()) == fProtonCode ||
                TMath::Abs(pdg->PdgCode()) == fNeutronCode ||
                TMath::Abs(pdg->PdgCode()) == fLambdaCode ||
                TMath::Abs(pdg->PdgCode()) == fXiCode ||
                TMath::Abs(pdg->PdgCode()) == fXi0Code ||
                TMath::Abs(pdg->PdgCode()) == fOmegaCode
	       )
            {
	      if (pdg->PdgCode() > 0) { particleMassPart = - protonMass;}
	      if (pdg->PdgCode() < 0) { particleMassPart = protonMass;}
	    }
	    Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
	    	      
	    if (pdg->PdgCode() == fProtonCode || pdg->PdgCode() == fAntiProtonCode)
	      {
		fProtonEt += et;
	      }
	    if (pdg->PdgCode() == fPiPlusCode || pdg->PdgCode() == fPiMinusCode)
	      {
		fPionEt += et;
	      }
	    if (pdg->PdgCode() == fKPlusCode || pdg->PdgCode() == fKMinusCode)
	      {
		fChargedKaonEt += et;
	      }
	    if (pdg->PdgCode() == fMuPlusCode || pdg->PdgCode() == fMuMinusCode)
	      {
		fMuonEt += et;
	      }
	    if (pdg->PdgCode() == fEPlusCode || pdg->PdgCode() == fEMinusCode)
{'fraction_non_alphanumeric': 0.142996965841245, 'fraction_numerical': 0.003915043554859548, 'mean_word_length': 3.305941845764855, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 129, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3208808', 'n_tokens_mistral': 3628, 'n_tokens_neox': 3358, 'n_words': 681}