Double check if SM is running added. Some redundant output removed from SM
[u/mrichter/AliRoot.git] / MFT / AliAnalysisTaskMFTExample.cxx
CommitLineData
d30a0b7b 1#include "TH1D.h"
2#include "AliAODEvent.h"
3#include "AliAODMCParticle.h"
4#include "AliAODMCHeader.h"
5#include "AliAnalysisTaskMFTExample.h"
6#include "TDatabasePDG.h"
7#include "AliMUONTrackParam.h"
8#include "AliMUONTrackExtrap.h"
9#include "AliAODDimuon.h"
10#include "AliAODTrack.h"
11#include "AliAODHeader.h"
12#include "TClonesArray.h"
13#include "AliMFTConstants.h"
bffc7f8c 14#include "AliMFTAnalysisTools.h"
d30a0b7b 15#include "TRandom.h"
2cf4030c 16#include "TList.h"
d30a0b7b 17
18ClassImp(AliAnalysisTaskMFTExample)
19
20//====================================================================================================================================================
21
22AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() :
23 AliAnalysisTaskSE(),
24 fVertexMode(0),
2cf4030c 25 fHistogramList(0),
d30a0b7b 26 fHistPtSingleMuons(0),
27 fHistPtSingleMuonsFromJpsi(0),
28 fHistPtDimuonsOS(0),
29 fHistMassDimuonsOS(0),
30 fHistPtDimuonsJpsi(0),
31 fHistMassDimuonsJpsi(0),
32 fHistResidualXVtxJpsi(0),
33 fHistResidualYVtxJpsi(0),
34 fHistResidualZVtxJpsi(0) {
35
36 // Default constructor
37
38 for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
39 fVtxResolutionITS[0] = 5.e-4;
40 fVtxResolutionITS[1] = 5.e-4;
41 fVtxResolutionITS[2] = 4.e-4;
42
43}
44
45//====================================================================================================================================================
46
47AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) :
48 AliAnalysisTaskSE(name),
49 fVertexMode(0),
2cf4030c 50 fHistogramList(0),
d30a0b7b 51 fHistPtSingleMuons(0),
52 fHistPtSingleMuonsFromJpsi(0),
53 fHistPtDimuonsOS(0),
54 fHistMassDimuonsOS(0),
55 fHistPtDimuonsJpsi(0),
56 fHistMassDimuonsJpsi(0),
57 fHistResidualXVtxJpsi(0),
58 fHistResidualYVtxJpsi(0),
59 fHistResidualZVtxJpsi(0) {
60
61 // Constructor
62
63 for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
64 fVtxResolutionITS[0] = 5.e-4;
65 fVtxResolutionITS[1] = 5.e-4;
66 fVtxResolutionITS[2] = 4.e-4;
67
68 // Define input and output slots here
2cf4030c 69 DefineOutput(1, TList::Class());
d30a0b7b 70
71}
72
73//====================================================================================================================================================
74
75void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {
76
77 // Called once
78
2cf4030c 79 fHistogramList = new TList();
80 fHistogramList->SetOwner(kTRUE);
81
d30a0b7b 82 fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
83 fHistPtSingleMuons -> SetXTitle("p_{T} [GeV/c]");
84 fHistPtSingleMuons -> Sumw2();
85
86 fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
87 fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T} [GeV/c]");
88 fHistPtSingleMuonsFromJpsi -> Sumw2();
89
90 fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
91 fHistMassDimuonsOS -> SetXTitle("Mass [GeV/c^{2}]");
92 fHistMassDimuonsOS -> Sumw2();
93
94 fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
95 fHistMassDimuonsJpsi -> SetXTitle("Mass [GeV/c^{2}]");
96 fHistMassDimuonsJpsi -> Sumw2();
97
98 fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
99 fHistPtDimuonsOS -> SetXTitle("p_{T} [GeV/c]");
100 fHistPtDimuonsOS -> Sumw2();
101
102 fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
103 fHistPtDimuonsJpsi -> SetXTitle("p_{T} [GeV/c]");
104 fHistPtDimuonsJpsi -> Sumw2();
105
106 fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
107 fHistResidualXVtxJpsi -> SetXTitle("Residual [#mum]");
108 fHistResidualXVtxJpsi -> Sumw2();
109
110 fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
111 fHistResidualYVtxJpsi -> SetXTitle("Residual [#mum]");
112 fHistResidualYVtxJpsi -> Sumw2();
113
114 fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
115 fHistResidualZVtxJpsi -> SetXTitle("Residual [#mum]");
116 fHistResidualZVtxJpsi -> Sumw2();
117
2cf4030c 118 fHistogramList->Add(fHistPtSingleMuons);
119 fHistogramList->Add(fHistPtSingleMuonsFromJpsi);
120 fHistogramList->Add(fHistMassDimuonsOS);
121 fHistogramList->Add(fHistMassDimuonsJpsi);
122 fHistogramList->Add(fHistPtDimuonsOS);
123 fHistogramList->Add(fHistPtDimuonsJpsi);
124 fHistogramList->Add(fHistResidualXVtxJpsi);
125 fHistogramList->Add(fHistResidualYVtxJpsi);
126 fHistogramList->Add(fHistResidualZVtxJpsi);
127
128 PostData(1, fHistogramList);
129
d30a0b7b 130}
131
132//====================================================================================================================================================
133
134void AliAnalysisTaskMFTExample::UserExec(Option_t *) {
135
136 // Main loop
137 // Called for each event
138
139 AliAODEvent *aodEv = dynamic_cast<AliAODEvent*>(InputEvent());
140 if (!aodEv) return;
141
142 aodEv->GetHeader()->InitMagneticField();
143 AliMUONTrackExtrap::SetField();
144
145 TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
146
147 // Getting primary vertex, either from the generation or from the reconstruction -------------------
148
149 if (fVertexMode == kGenerated) {
150 AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
151 mcHeader->GetVertex(fPrimaryVertex);
152 for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]);
153 }
154 else if (fVertexMode == kReconstructed) {
155 aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex);
156 }
157
158 // -------------------------------------------------------------------------------------------------
159
160 AliAODTrack *recMuon1=0, *recMuon2=0;
161 AliAODMCParticle *mcMuon1=0, *mcMuon2=0;
162 AliAODMCParticle *mcMother1=0;
163 Bool_t isMuon1FromJpsi=0;
164
165 // Loop over MUON+MFT muons
166
167 //--------------------------------------------------------------------------------------------------
168
169 for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) {
170
171 if (!(aodEv->GetTrack(iTrack)->IsMuonGlobalTrack())) continue;
172
173 recMuon1 = aodEv->GetTrack(iTrack);
174 isMuon1FromJpsi = kFALSE;
175
176 // pt all muons
177 fHistPtSingleMuons -> Fill(recMuon1->Pt());
178
179 // pt muons from J/psi
180 if (recMuon1->GetLabel() >= 0) {
181 mcMuon1 = (AliAODMCParticle*) stackMC->At(recMuon1->GetLabel());
182 if (mcMuon1) {
183 if (mcMuon1->GetMother()>=0) {
184 mcMother1 = (AliAODMCParticle*) stackMC->At(mcMuon1->GetMother());
185 if (mcMother1->PdgCode()==443) {
186 isMuon1FromJpsi = kTRUE;
187 fHistPtSingleMuonsFromJpsi -> Fill(recMuon1->Pt());
188 }
189 }
190 }
191 }
192
193 for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
194
195 if (!(aodEv->GetTrack(jTrack)->IsMuonGlobalTrack())) continue;
196
197 recMuon2 = aodEv->GetTrack(jTrack);
198
199 AliAODDimuon *dimuon = new AliAODDimuon(recMuon1, recMuon2);
0ebb458a 200 if (dimuon->Charge()) {
201 delete dimuon;
202 continue;
203 }
d30a0b7b 204
205 // pt and mass all OS dimuons
206 fHistPtDimuonsOS -> Fill(dimuon->Pt());
207 fHistMassDimuonsOS -> Fill(dimuon->Mass());
208
0ebb458a 209 delete dimuon;
210
d30a0b7b 211 // pt and mass J/psi dimuons
212 if (!isMuon1FromJpsi) continue;
213 if (recMuon2->GetLabel() >= 0) {
214 mcMuon2 = (AliAODMCParticle*) stackMC->At(recMuon2->GetLabel());
215 if (mcMuon2) {
216 if (mcMuon2->GetMother() == mcMuon1->GetMother()) {
0ebb458a 217 AliAODDimuon *dimuonJpsi = new AliAODDimuon;
218 dimuonJpsi->SetMuons(recMuon1,recMuon2);
bffc7f8c 219 Double_t pca[3]={0};
220 Double_t pcaQuality=0;
d30a0b7b 221 TLorentzVector kinem(0,0,0,0);
0ebb458a 222 if (!AliMFTAnalysisTools::CalculatePCA(dimuonJpsi, pca, pcaQuality, kinem)) {
223 delete dimuonJpsi;
224 continue;
225 }
d30a0b7b 226 fHistPtDimuonsJpsi -> Fill(kinem.Pt());
227 fHistMassDimuonsJpsi -> Fill(kinem.M());
bffc7f8c 228 fHistResidualXVtxJpsi -> Fill(1.e4*(pca[0] - fPrimaryVertex[0]));
229 fHistResidualYVtxJpsi -> Fill(1.e4*(pca[1] - fPrimaryVertex[1]));
230 fHistResidualZVtxJpsi -> Fill(1.e4*(pca[2] - fPrimaryVertex[2]));
0ebb458a 231 delete dimuonJpsi;
d30a0b7b 232 }
233 }
234 }
235
236 } // end of loop on 2nd muon
237
238 } // end of loop on 1st muon
239
240 //--------------------------------------------------------------------------------------------------
241
2cf4030c 242 PostData(1, fHistogramList);
d30a0b7b 243
244}
245
246//====================================================================================================================================================
247
248void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
249
250 // Draw result to the screen
251 // Called once at the end of the query
252
253}
254
255//====================================================================================================================================================