PWG4/totEt/AliAnalysisEtMonteCarlo.cxx

   1 #include "AliAnalysisEtMonteCarlo.h"
   2 #include "AliAnalysisEtCuts.h"
   3
   4 #include "AliStack.h"
   5 #include "AliMCEvent.h"
   6
   7 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
   8 {
   9      ResetEventValues();
  10
  11     // Get us an mc event
  12     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
  13
  14     // Let's play with the stack!
  15     AliStack *stack = event->Stack();
  16
  17     Int_t nPrim = stack->GetNtrack();
  18
  19     Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
  20
  21     for (Int_t iPart = 0; iPart < nPrim; iPart++)
  22     {
  23
  24         TParticle *part = stack->Particle(iPart);
  25
  26         if (!part)
  27         {
  28             Printf("ERROR: Could not get particle %d", iPart);
  29             continue;
  30         }
  31
  32         TParticlePDG *pc = part->GetPDG(0);
  33
  34         // Check if it is a primary particle
  35         if (!stack->IsPhysicalPrimary(iPart)) continue;
  36
  37         // Check for reasonable (for now neutral and singly charged) charge on the particle
  38         //TODO:Maybe not only singly charged?
  39         if (TMath::Abs(pc->Charge()) != EtMonteCarloCuts::kSingleChargedParticle && pc->Charge() != EtMonteCarloCuts::kNeutralParticle) continue;
  40
  41         fMultiplicity++;
  42
  43         if (TMath::Abs(part->Eta()) < fEtaCut)
  44         {
  45
  46             TParticlePDG *pdgCode =  part->GetPDG(0);
  47             if (
  48                 TMath::Abs(pdgCode->PdgCode()) == fPdgDB->GetParticle("proton")->PdgCode() ||
  49                 TMath::Abs(pdgCode->PdgCode()) == fPdgDB->GetParticle("neutron")->PdgCode() ||
  50                 TMath::Abs(pdgCode->PdgCode()) == fPdgDB->GetParticle("Lambda0")->PdgCode() ||
  51                 TMath::Abs(pdgCode->PdgCode()) == fPdgDB->GetParticle("Xi-")->PdgCode() ||
  52                 TMath::Abs(pdgCode->PdgCode()) == fPdgDB->GetParticle("Xi0")->PdgCode() ||
  53                 TMath::Abs(pdgCode->PdgCode()) == fPdgDB->GetParticle("Omega-")->PdgCode()
  54                )
  55             {
  56                 particleMassPart = -TMath::Sign(pdgCode->PdgCode(), pdgCode->PdgCode())*pdgCode->Mass();
  57             }
  58
  59             if (pdgCode->Charge() == EtMonteCarloCuts::kNeutralParticle)
  60             {
  61                fNeutralMultiplicity++;
  62                 fTotNeutralEt += part->Energy()*TMath::Sin(part->Theta());
  63
  64                 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
  65                 {
  66                     fTotNeutralEtAcc += part->Energy()*TMath::Sin(part->Theta());
  67                     fTotEtAcc += part->Energy()*TMath::Sin(part->Theta());
  68                 }
  69             }
  70             else if (pdgCode->Charge() != EtMonteCarloCuts::kNeutralParticle)
  71             {
  72                fChargedMultiplicity++;
  73                 fTotChargedEt += part->Energy()*TMath::Sin(part->Theta());
  74                 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
  75                 {
  76                     fTotChargedEtAcc += part->Energy()*TMath::Sin(part->Theta());
  77                     fTotEtAcc += part->Energy()*TMath::Sin(part->Theta());
  78                 }
  79             }
  80         }
  81     }
  82
  83     fTotNeutralEtAcc = fTotNeutralEt;
  84     fTotEt = fTotChargedEt + fTotNeutralEt;
  85     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
  86
  87     FillHistograms();
  88
  89     return 0;
  90 }
  91
  92 void AliAnalysisEtMonteCarlo::Init()
  93 {
  94
  95     AliAnalysisEt::Init();
  96
  97     fVertexXCut = EtReconstructedCuts::kVertexXCut;
  98     fVertexYCut = EtReconstructedCuts::kVertexYCut;
  99     fVertexZCut = EtReconstructedCuts::kVertexZCut;
 100     fIPxyCut = EtReconstructedCuts::kIPxyCut;
 101     fIPzCut = EtReconstructedCuts::kIPzCut;
 102
 103 }