PWG4/totEt/AliAnalysisEtMonteCarlo.cxx

   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
  10 #include "AliAnalysisEtMonteCarlo.h"
  11 #include "AliAnalysisEtCuts.h"
  12 #include "AliESDtrack.h"
  13 #include "AliStack.h"
  14 #include "AliMCEvent.h"
  15 #include "TH2F.h"
  16 #include "TParticle.h"
  17
  18 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
  19 { // analyse MC event
  20      ResetEventValues();
  21
  22     // Get us an mc event
  23     AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
  24
  25     // Let's play with the stack!
  26     AliStack *stack = event->Stack();
  27
  28     Int_t nPrim = stack->GetNtrack();
  29
  30     Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
  31
  32     for (Int_t iPart = 0; iPart < nPrim; iPart++)
  33     {
  34
  35         TParticle *part = stack->Particle(iPart);
  36
  37         if (!part)
  38         {
  39             Printf("ERROR: Could not get particle %d", iPart);
  40             continue;
  41         }
  42
  43         TParticlePDG *pc = part->GetPDG(0);
  44
  45         // Check if it is a primary particle
  46         if (!stack->IsPhysicalPrimary(iPart)) continue;
  47
  48         // 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
  49
  50         // Check for reasonable (for now neutral and singly charged) charge on the particle
  51         //TODO:Maybe not only singly charged?
  52         if (TMath::Abs(pc->Charge()) != EtMonteCarloCuts::kSingleChargedParticle && pc->Charge() != EtMonteCarloCuts::kNeutralParticle) continue;
  53
  54         fMultiplicity++;
  55
  56         if (TMath::Abs(part->Eta()) < fEtaCut)
  57         {
  58
  59           TParticlePDG *pdgCode =  part->GetPDG(0);
  60             if (
  61                 TMath::Abs(pdgCode->PdgCode()) == fProtonCode ||
  62                 TMath::Abs(pdgCode->PdgCode()) == fNeutronCode ||
  63                 TMath::Abs(pdgCode->PdgCode()) == fLambdaCode ||
  64                 TMath::Abs(pdgCode->PdgCode()) == fXiCode ||
  65                 TMath::Abs(pdgCode->PdgCode()) == fXi0Code ||
  66                 TMath::Abs(pdgCode->PdgCode()) == fOmegaCode
  67                )
  68             {
  69                 particleMassPart = -TMath::Sign(pdgCode->PdgCode(), pdgCode->PdgCode())*pdgCode->Mass();
  70             }
  71
  72             if (pdgCode->Charge() == EtMonteCarloCuts::kNeutralParticle)
  73             {
  74                fNeutralMultiplicity++;
  75                 fTotNeutralEt += part->Energy()*TMath::Sin(part->Theta());
  76
  77                 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
  78                 {
  79                     fTotNeutralEtAcc += part->Energy()*TMath::Sin(part->Theta());
  80                     fTotEtAcc += part->Energy()*TMath::Sin(part->Theta());
  81                 }
  82             }
  83             else if (pdgCode->Charge() != EtMonteCarloCuts::kNeutralParticle)
  84             {
  85                fChargedMultiplicity++;
  86                 fTotChargedEt += part->Energy()*TMath::Sin(part->Theta());
  87                 if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
  88                 {
  89                     fTotChargedEtAcc += part->Energy()*TMath::Sin(part->Theta());
  90                     fTotEtAcc += part->Energy()*TMath::Sin(part->Theta());
  91                 }
  92
  93                 //        if (TrackHitsCalorimeter(part, event->GetMagneticField()))
  94                 if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
  95                   {
  96                     if (pdgCode->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
  97                     else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
  98                   }
  99             }
 100         }
 101     }
 102
 103     fTotNeutralEtAcc = fTotNeutralEt;
 104     fTotEt = fTotChargedEt + fTotNeutralEt;
 105     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
 106
 107     FillHistograms();
 108
 109     return 0;
 110 }
 111
 112 void AliAnalysisEtMonteCarlo::Init()
 113 {
 114     AliAnalysisEt::Init();
 115 }
 116
 117 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
 118 {
 119   //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
 120    AliESDtrack *esdTrack = new AliESDtrack(part);
 121    // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 122
 123     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 124
 125     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 126
 127     bool status = prop &&
 128       TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
 129       esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
 130       esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
 131     delete esdTrack;
 132
 133     return status;
 134 }
 135