[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"
#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
     if(!ev){
            Printf("ERROR: Event does not exist");   
	    return 0;
     }
     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);

    Double_t protonMass =fgProtonMass;

    // Hijing header
    AliGenEventHeader* genHeader = event->GenEventHeader();
    if(!genHeader){
            Printf("ERROR: Event generation header does not exist");   
	    return 0;
    }
    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()) == fgProtonCode ||
                TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
                TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
                TMath::Abs(pdg->PdgCode()) == fgXiCode ||
                TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
                TMath::Abs(pdg->PdgCode()) == fgOmegaCode
	       )
            {
	      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() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
	      {
		fProtonEt += et;
	      }
	    if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
	      {
		fPionEt += et;
	      }
	    if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
	      {
		fChargedKaonEt += et;
	      }
	    if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
	      {
		fMuonEt += et;
	      }
	    if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
	      {
		fElectronEt += et;
	      }

	    // some neutrals also
	    if(pdg->PdgCode() == fgNeutronCode)
	    {
	      fNeutronEt += et;
	    }
            if(pdg->PdgCode() == fgAntiNeutronCode)
	    {
	      fAntiNeutronEt += et;
	    }
	    if(pdg->PdgCode() == fgGammaCode)
	    {
	      fGammaEt += et;
	    }

            if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
            {
	       fNeutralMultiplicity++;
	       fTotNeutralEt += et;

                if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
                {
		  fTotNeutralEtAcc += et;
		  fTotEtAcc += et;
                }
            }
            else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
            {
	       fChargedMultiplicity++;
	       fTotChargedEt += et;
                if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
                {
		  fTotChargedEtAcc += et;
		  fTotEtAcc += et;


		    if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
		      {
			fProtonEtAcc += et;
		      }
		    if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
		      {
			fPionEtAcc += et;
		      }
		    if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
		      {
			fChargedKaonEtAcc += et;
		      }
		    if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
		      {
			fMuonEtAcc += et;
		      }
		    if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
		      {
			fElectronEtAcc += et;
		      }
		    
                }

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

    return 0;    
}

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

void AliAnalysisEtMonteCarlo::ResetEventValues()
{ // reset event values
  AliAnalysisEt::ResetEventValues();

  // collision geometry defaults for p+p:
  fImpactParameter = 0;
  fNcoll = 1;
  fNpart = 2;  
}

void AliAnalysisEtMonteCarlo::CreateHistograms()
{ // histogram related additions
  AliAnalysisEt::CreateHistograms();
  if (fTree) {
    fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
    fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
    fTree->Branch("fNpart",&fNpart,"fNpart/I");
  }
}

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;
}
Commit	Line	Data
cf6522d1	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Base class for MC analysis
	4	// - MC output
	5	// implementation file
	6	//
	7	//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
	8	//_________________________________________________________________________
	9
2fbf38ac	10	#include "AliAnalysisEtMonteCarlo.h"
2fbf38ac	11	#include "AliAnalysisEtCuts.h"
13b0d3c1	12	#include "AliESDtrack.h"
2fbf38ac	13	#include "AliStack.h"
2fbf38ac	14	#include "AliMCEvent.h"
13b0d3c1	15	#include "TH2F.h"
cf6522d1	16	#include "TParticle.h"
ce546038	17	#include "AliGenHijingEventHeader.h"
	18	#include "AliGenPythiaEventHeader.h"
	19
16abb579	20	using namespace std;
	21
	22	ClassImp(AliAnalysisEtMonteCarlo);
	23
	24
ce546038	25	// ctor
	26	AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo() :
	27	AliAnalysisEt()
	28	,fImpactParameter(0)
	29	,fNcoll(0)
	30	,fNpart(0)
	31	{
	32	}
	33
	34	// dtor
	35	AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
	36	{
	37	}
2fbf38ac	38
2fbf38ac	39	Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
cf6522d1	40	{ // analyse MC event
2fbf38ac	41	ResetEventValues();
	42
	43	// Get us an mc event
6a0df78a	44	if(!ev){
	45	Printf("ERROR: Event does not exist");
	46	return 0;
	47	}
	48	AliMCEvent event = dynamic_cast<AliMCEvent>(ev);
2fbf38ac	49
7d2d1773	50	Double_t protonMass =fgProtonMass;
91a79ea6	51
ce546038	52	// Hijing header
ce546038	53	AliGenEventHeader* genHeader = event->GenEventHeader();
f4675a17	54	if(!genHeader){
	55	Printf("ERROR: Event generation header does not exist");
	56	return 0;
	57	}
ce546038	58	AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
	59	if (hijingGenHeader) {
	60	fImpactParameter = hijingGenHeader->ImpactParameter();
	61	fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
	62	fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
	63	/*
	64	printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
	65	hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
	66	printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
	67	hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
	68	printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
	69	hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
	70	printf("NN %d NNw %d NwN %d, NwNw %d\n",
	71	hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
	72	*/
	73	}
	74
	75	/* // placeholder if we want to get some Pythia info later
	76	AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
	77	if (pythiaGenHeader) { // not Hijing; try with Pythia
	78	printf("Pythia: ProcessType %d GetPtHard %g \n",
	79	pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
	80	}
	81	*/
	82
2fbf38ac	83	// Let's play with the stack!
	84	AliStack *stack = event->Stack();
	85
	86	Int_t nPrim = stack->GetNtrack();
	87
2fbf38ac	88	for (Int_t iPart = 0; iPart < nPrim; iPart++)
	89	{
	90
	91	TParticle *part = stack->Particle(iPart);
	92
	93	if (!part)
	94	{
	95	Printf("ERROR: Could not get particle %d", iPart);
	96	continue;
	97	}
	98
7b873813	99	TParticlePDG *pdg = part->GetPDG(0);
	100
	101	Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
2fbf38ac	102
	103	// Check if it is a primary particle
	104	if (!stack->IsPhysicalPrimary(iPart)) continue;
	105
7b873813	106	// 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
cf6522d1	107
2fbf38ac	108	// Check for reasonable (for now neutral and singly charged) charge on the particle
2fbf38ac	109	//TODO:Maybe not only singly charged?
464aa50c	110	if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
2fbf38ac	111
	112	fMultiplicity++;
	113
4998becf	114	if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
2fbf38ac	115	{
2fbf38ac	116
2fbf38ac	117	if (
7d2d1773	118	TMath::Abs(pdg->PdgCode()) == fgProtonCode \|\|
	119	TMath::Abs(pdg->PdgCode()) == fgNeutronCode \|\|
	120	TMath::Abs(pdg->PdgCode()) == fgLambdaCode \|\|
	121	TMath::Abs(pdg->PdgCode()) == fgXiCode \|\|
	122	TMath::Abs(pdg->PdgCode()) == fgXi0Code \|\|
	123	TMath::Abs(pdg->PdgCode()) == fgOmegaCode
2fbf38ac	124	)
2fbf38ac	125	{
91a79ea6	126	if (pdg->PdgCode() > 0) { particleMassPart = - protonMass;}
91a79ea6	127	if (pdg->PdgCode() < 0) { particleMassPart = protonMass;}
43056f1b	128	}
7b873813	129	Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
7b873813	130
7d2d1773	131	if (pdg->PdgCode() == fgProtonCode \|\| pdg->PdgCode() == fgAntiProtonCode)
7b873813	132	{
	133	fProtonEt += et;
	134	}
7d2d1773	135	if (pdg->PdgCode() == fgPiPlusCode \|\| pdg->PdgCode() == fgPiMinusCode)
7b873813	136	{
	137	fPionEt += et;
	138	}
7d2d1773	139	if (pdg->PdgCode() == fgKPlusCode \|\| pdg->PdgCode() == fgKMinusCode)
7b873813	140	{
	141	fChargedKaonEt += et;
	142	}
7d2d1773	143	if (pdg->PdgCode() == fgMuPlusCode \|\| pdg->PdgCode() == fgMuMinusCode)
7b873813	144	{
	145	fMuonEt += et;
	146	}
7d2d1773	147	if (pdg->PdgCode() == fgEPlusCode \|\| pdg->PdgCode() == fgEMinusCode)
7b873813	148	{
	149	fElectronEt += et;
	150	}
	151
	152	// some neutrals also
7d2d1773	153	if(pdg->PdgCode() == fgNeutronCode)
43056f1b	154	{
7b873813	155	fNeutronEt += et;
43056f1b	156	}
7d2d1773	157	if(pdg->PdgCode() == fgAntiNeutronCode)
43056f1b	158	{
7b873813	159	fAntiNeutronEt += et;
43056f1b	160	}
7d2d1773	161	if(pdg->PdgCode() == fgGammaCode)
43056f1b	162	{
7b873813	163	fGammaEt += et;
43056f1b	164	}
7b873813	165
464aa50c	166	if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
2fbf38ac	167	{
2fbf38ac	168	fNeutralMultiplicity++;
7b873813	169	fTotNeutralEt += et;
2fbf38ac	170
	171	if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
	172	{
7b873813	173	fTotNeutralEtAcc += et;
7b873813	174	fTotEtAcc += et;
2fbf38ac	175	}
2fbf38ac	176	}
b527722d	177	else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
2fbf38ac	178	{
2fbf38ac	179	fChargedMultiplicity++;
7b873813	180	fTotChargedEt += et;
2fbf38ac	181	if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
2fbf38ac	182	{
7b873813	183	fTotChargedEtAcc += et;
	184	fTotEtAcc += et;
	185
	186
7d2d1773	187	if (pdg->PdgCode() == fgProtonCode \|\| pdg->PdgCode() == fgAntiProtonCode)
7b873813	188	{
	189	fProtonEtAcc += et;
	190	}
7d2d1773	191	if (pdg->PdgCode() == fgPiPlusCode \|\| pdg->PdgCode() == fgPiMinusCode)
7b873813	192	{
	193	fPionEtAcc += et;
	194	}
7d2d1773	195	if (pdg->PdgCode() == fgKPlusCode \|\| pdg->PdgCode() == fgKMinusCode)
7b873813	196	{
	197	fChargedKaonEtAcc += et;
	198	}
7d2d1773	199	if (pdg->PdgCode() == fgMuPlusCode \|\| pdg->PdgCode() == fgMuMinusCode)
7b873813	200	{
	201	fMuonEtAcc += et;
	202	}
7d2d1773	203	if (pdg->PdgCode() == fgEPlusCode \|\| pdg->PdgCode() == fgEMinusCode)
7b873813	204	{
	205	fElectronEtAcc += et;
	206	}
	207
2fbf38ac	208	}
13b0d3c1	209
	210	// if (TrackHitsCalorimeter(part, event->GetMagneticField()))
	211	if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
	212	{
7b873813	213	if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
13b0d3c1	214	else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
	215	}
	216	}
	217	}
2fbf38ac	218	}
2fbf38ac	219
2fbf38ac	220	fTotEt = fTotChargedEt + fTotNeutralEt;
	221	fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
	222
	223	FillHistograms();
	224
	225	return 0;
	226	}
	227
	228	void AliAnalysisEtMonteCarlo::Init()
ce546038	229	{ // init
2fbf38ac	230	AliAnalysisEt::Init();
2fbf38ac	231	}
13b0d3c1	232
ce546038	233	void AliAnalysisEtMonteCarlo::ResetEventValues()
	234	{ // reset event values
	235	AliAnalysisEt::ResetEventValues();
	236
	237	// collision geometry defaults for p+p:
	238	fImpactParameter = 0;
	239	fNcoll = 1;
	240	fNpart = 2;
	241	}
	242
	243	void AliAnalysisEtMonteCarlo::CreateHistograms()
	244	{ // histogram related additions
	245	AliAnalysisEt::CreateHistograms();
	246	if (fTree) {
	247	fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
	248	fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
	249	fTree->Branch("fNpart",&fNpart,"fNpart/I");
	250	}
	251	}
	252
13b0d3c1	253	bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
	254	{
	255	// printf(" TrackHitsCalorimeter - magField %f\n", magField);
	256	AliESDtrack *esdTrack = new AliESDtrack(part);
	257	// Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
	258
	259	Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
	260
	261	// if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
	262
	263	bool status = prop &&
	264	TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
	265	esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
	266	esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
	267	delete esdTrack;
	268
	269	return status;
	270	}
	271