PWG4/totEt/AliAnalysisEtMonteCarlo.cxx

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