]>
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" |
11 | #include "AliAnalysisEtCuts.h" | |
13b0d3c1 | 12 | #include "AliESDtrack.h" |
2fbf38ac | 13 | #include "AliStack.h" |
14 | #include "AliMCEvent.h" | |
13b0d3c1 | 15 | #include "TH2F.h" |
cf6522d1 | 16 | #include "TParticle.h" |
2fbf38ac | 17 | |
18 | Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev) | |
cf6522d1 | 19 | { // analyse MC event |
2fbf38ac | 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 | ||
cf6522d1 | 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 | ||
2fbf38ac | 50 | // Check for reasonable (for now neutral and singly charged) charge on the particle |
51 | //TODO:Maybe not only singly charged? | |
4998becf | 52 | if (TMath::Abs(pc->Charge()) != fCuts->GetMonteCarloSingleChargedParticle() && pc->Charge() != fCuts->GetMonteCarloNeutralParticle()) continue; |
2fbf38ac | 53 | |
54 | fMultiplicity++; | |
55 | ||
4998becf | 56 | if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut()) |
2fbf38ac | 57 | { |
58 | ||
13b0d3c1 | 59 | TParticlePDG *pdgCode = part->GetPDG(0); |
2fbf38ac | 60 | if ( |
8985c506 | 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 | |
2fbf38ac | 67 | ) |
68 | { | |
69 | particleMassPart = -TMath::Sign(pdgCode->PdgCode(), pdgCode->PdgCode())*pdgCode->Mass(); | |
70 | } | |
71 | ||
4998becf | 72 | if (pdgCode->Charge() == fCuts->GetMonteCarloNeutralParticle() ) |
2fbf38ac | 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 | } | |
4998becf | 83 | else if (pdgCode->Charge() != fCuts->GetMonteCarloNeutralParticle() ) |
2fbf38ac | 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 | } | |
13b0d3c1 | 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 | } | |
2fbf38ac | 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 | { | |
2fbf38ac | 114 | AliAnalysisEt::Init(); |
2fbf38ac | 115 | } |
13b0d3c1 | 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 |