totEt cuts into their own proper class
[u/mrichter/AliRoot.git] / 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()) != fCuts->GetMonteCarloSingleChargedParticle() && pc->Charge() != fCuts->GetMonteCarloNeutralParticle()) continue;
53
54         fMultiplicity++;
55
56         if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
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() == fCuts->GetMonteCarloNeutralParticle() )
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() != fCuts->GetMonteCarloNeutralParticle() )
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